GIT: relocate treeDataEdge, treeDataPoint into OpenFOAM

- similar to treeDataCell, doesn't need any meshTools components
- AABBTree under algorithms (like indexedOctree)

src/OpenFOAM/algorithms/AABBTree/AABBTree.C

@@ -0,0 +1,539 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2015 OpenFOAM Foundation
+    Copyright (C) 2016-2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "AABBTree.H"
+#include "bitSet.H"
+
+template<class Type>
+Foam::scalar Foam::AABBTree<Type>::tolerance_ = 1e-4;
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class Type>
+void Foam::AABBTree<Type>::writeOBJ
+(
+    const bool writeLinesOnly,
+    const treeBoundBox& bb,
+    label& vertI,
+    Ostream& os
+) const
+{
+    const pointField pts(bb.points());
+
+    for (const point& p : pts)
+    {
+        os << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
+    }
+
+    if (writeLinesOnly)
+    {
+        for (const edge& e : bb.edges)
+        {
+            os  << "l " << e[0] + vertI + 1 << ' ' << e[1] + vertI + 1 << nl;
+        }
+    }
+    else
+    {
+        for (const face& f : bb.faces)
+        {
+            os  << 'f';
+            for (const label fpi : f)
+            {
+                os  << ' ' << fpi + vertI + 1;
+            }
+            os  << nl;
+        }
+    }
+
+    vertI += pts.size();
+}
+
+
+template<class Type>
+void Foam::AABBTree<Type>::writeOBJ
+(
+    const bool leavesOnly,
+    const bool writeLinesOnly,
+    const treeBoundBox& bb,
+    const label nodeI,
+    const List<Pair<treeBoundBox>>& bbs,
+    const List<Pair<label>>& nodes,
+    label& vertI,
+    Ostream& os
+) const
+{
+    if (!leavesOnly || nodeI < 0)
+    {
+        writeOBJ(writeLinesOnly, bb, vertI, os);
+    }
+
+    // recurse to find leaves
+    if (nodeI >= 0)
+    {
+        writeOBJ
+        (
+            leavesOnly,
+            writeLinesOnly,
+            bbs[nodeI].first(),
+            nodes[nodeI].first(),
+            bbs,
+            nodes,
+            vertI,
+            os
+        );
+        writeOBJ
+        (
+            leavesOnly,
+            writeLinesOnly,
+            bbs[nodeI].second(),
+            nodes[nodeI].second(),
+            bbs,
+            nodes,
+            vertI,
+            os
+        );
+    }
+}
+
+
+template<class Type>
+void Foam::AABBTree<Type>::createBoxes
+(
+    const bool equalBinSize,
+    const label level,
+    const List<Type>& objects,
+    const pointField& points,
+    const DynamicList<label>& objectIDs,
+    const treeBoundBox& bb,
+    const label nodeI,
+
+    DynamicList<Pair<treeBoundBox>>& bbs,
+    DynamicList<labelPair>& nodes,
+    DynamicList<labelList>& addressing
+) const
+{
+    const vector span = bb.span();
+
+    // Determine which direction to divide the box
+
+    direction maxDir = 0;
+    scalar maxSpan = span[maxDir];
+    for (label dirI = 1; dirI < 3; ++dirI)
+    {
+        if (span[dirI] > maxSpan)
+        {
+            maxSpan = span[dirI];
+            maxDir = dirI;
+        }
+    }
+
+
+    scalar divide;
+
+    if (equalBinSize)
+    {
+        // Pick up points used by this set of objects
+
+        bitSet isUsedPoint(points.size());
+        DynamicList<scalar> component(points.size());
+
+        for (const label objI : objectIDs)
+        {
+            const Type& obj = objects[objI];
+
+            for (const label pointI : obj)
+            {
+                if (isUsedPoint.set(pointI))
+                {
+                    component.append(points[pointI][maxDir]);
+                }
+            }
+        }
+
+        // Determine the median
+
+        Foam::sort(component);
+
+        divide = component[component.size()/2];
+    }
+    else
+    {
+        // Geometric middle
+        divide = bb.min()[maxDir] + 0.5*maxSpan;
+    }
+
+
+    scalar divMin = divide + tolerance_*maxSpan;
+    scalar divMax = divide - tolerance_*maxSpan;
+
+
+    // Assign the objects to min or max bin
+
+    DynamicList<label> minBinObjectIDs(objectIDs.size());
+    treeBoundBox minBb(boundBox::invertedBox);
+
+    DynamicList<label> maxBinObjectIDs(objectIDs.size());
+    treeBoundBox maxBb(boundBox::invertedBox);
+
+    for (const label objI : objectIDs)
+    {
+        const Type& obj = objects[objI];
+
+        bool intoMin = false;
+        bool intoMax = false;
+
+        for (const label pointI : obj)
+        {
+            const point& pt = points[pointI];
+            if (pt[maxDir] < divMin)
+            {
+                intoMin = true;
+            }
+            if (pt[maxDir] > divMax)
+            {
+                intoMax = true;
+            }
+        }
+
+        // Note: object is inserted into both min/max child boxes (duplicated)
+        // if it crosses the bin boundaries
+        if (intoMin)
+        {
+            minBinObjectIDs.append(objI);
+            minBb.add(points, obj);
+        }
+        if (intoMax)
+        {
+            maxBinObjectIDs.append(objI);
+            maxBb.add(points, obj);
+        }
+    }
+
+    // Inflate box in case geometry reduces to 2-D
+    if (minBinObjectIDs.size())
+    {
+        minBb.inflate(0.01);
+    }
+    if (maxBinObjectIDs.size())
+    {
+        maxBb.inflate(0.01);
+    }
+
+    minBinObjectIDs.shrink();
+    maxBinObjectIDs.shrink();
+
+
+    label minI;
+    if (minBinObjectIDs.size() > minLeafSize_ && level < maxLevel_)
+    {
+        // New leaf
+        minI = nodes.size();
+        nodes.append(labelPair(-1, -1));
+    }
+    else
+    {
+        // Update existing leaf
+        minI = -addressing.size() - 1;
+        addressing.append(minBinObjectIDs);
+    }
+
+    label maxI;
+    if (maxBinObjectIDs.size() > minLeafSize_ && level < maxLevel_)
+    {
+        // New leaf
+        maxI = nodes.size();
+        nodes.append(labelPair(-1, -1));
+    }
+    else
+    {
+        // Update existing leaf
+        maxI = -addressing.size() - 1;
+        addressing.append(maxBinObjectIDs);
+    }
+
+    nodes(nodeI) = labelPair(minI, maxI);
+    bbs(nodeI) = Pair<treeBoundBox>(minBb, maxBb);
+
+    // Recurse
+    if (minI >= 0)
+    {
+        createBoxes
+        (
+            equalBinSize,
+            level + 1,
+            objects,
+            points,
+            minBinObjectIDs,
+            minBb,
+            minI,
+            bbs,
+            nodes,
+            addressing
+        );
+    }
+    if (maxI >= 0)
+    {
+        createBoxes
+        (
+            equalBinSize,
+            level + 1,
+            objects,
+            points,
+            maxBinObjectIDs,
+            maxBb,
+            maxI,
+            bbs,
+            nodes,
+            addressing
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::AABBTree<Type>::AABBTree()
+:
+    maxLevel_(0),
+    minLeafSize_(0),
+    boundBoxes_(),
+    addressing_()
+{}
+
+
+template<class Type>
+Foam::AABBTree<Type>::AABBTree
+(
+    const UList<Type>& objects,
+    const pointField& points,
+    const bool equalBinSize,
+    const label maxLevel,
+    const label minLeafSize
+)
+:
+    maxLevel_(maxLevel),
+    minLeafSize_(minLeafSize),
+    boundBoxes_(),
+    addressing_()
+{
+    if (objects.empty())
+    {
+        return;
+    }
+
+
+    DynamicList<Pair<treeBoundBox>> bbs(maxLevel);
+    DynamicList<labelPair> nodes(maxLevel);
+    DynamicList<labelList> addr(maxLevel);
+
+    nodes.append(labelPair(-1, -1));
+    treeBoundBox topBb(points);
+    topBb.inflate(0.01);
+
+    DynamicList<label> objectIDs(identity(objects.size()));
+
+    createBoxes
+    (
+        equalBinSize,
+        0,          // starting at top level
+        objects,
+        points,
+        objectIDs,
+        topBb,
+        0,          // starting node
+
+        bbs,
+        nodes,
+        addr
+    );
+
+
+    //{
+    //    OFstream os("tree.obj");
+    //    label vertI = 0;
+    //    writeOBJ
+    //    (
+    //        true,       // leavesOnly
+    //        false,      // writeLinesOnly
+    //
+    //        topBb,
+    //        0,
+    //        bbs,
+    //        nodes,
+    //        vertI,
+    //        os
+    //    );
+    //}
+
+
+    // transfer flattened tree to persistent storage
+    DynamicList<treeBoundBox> boundBoxes(2*bbs.size());
+    DynamicList<labelList> addressing(2*addr.size());
+
+    forAll(nodes, nodeI)
+    {
+        if (nodes[nodeI].first() < 0)
+        {
+            boundBoxes.append(bbs[nodeI].first());
+            addressing.append(addr[-(nodes[nodeI].first() + 1)]);
+        }
+        if (nodes[nodeI].second() < 0)
+        {
+            boundBoxes.append(bbs[nodeI].second());
+            addressing.append(addr[-(nodes[nodeI].second() + 1)]);
+        }
+    }
+
+    boundBoxes_.transfer(boundBoxes);
+    addressing_.transfer(addressing);
+
+
+    if (0)
+    {
+        bitSet checked(objects.size());
+        for (const auto& box : addressing_)
+        {
+            for (const auto& id : box)
+            {
+                checked.set(id);
+            }
+        }
+
+        const label unsetSize = checked.count(false);
+
+        if (unsetSize)
+        {
+            Info<< "*** Problem: IDs not set: " << unsetSize << endl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+const Foam::List<Foam::treeBoundBox>& Foam::AABBTree<Type>::boundBoxes() const
+{
+    return boundBoxes_;
+}
+
+
+template<class Type>
+const Foam::List<Foam::labelList>& Foam::AABBTree<Type>::addressing() const
+{
+    return addressing_;
+}
+
+
+template<class Type>
+bool Foam::AABBTree<Type>::pointInside(const point& pt) const
+{
+    for (const treeBoundBox& bb : boundBoxes_)
+    {
+        if (bb.contains(pt))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+template<class Type>
+bool Foam::AABBTree<Type>::overlaps(const boundBox& bbIn) const
+{
+    for (const treeBoundBox& bb : boundBoxes_)
+    {
+        if (bb.overlaps(bbIn))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Ostream& Foam::operator<<(Ostream& os, const AABBTree<Type>& tree)
+{
+    if (os.format() == IOstreamOption::ASCII)
+    {
+        os  << tree.maxLevel_ << token::SPACE
+            << tree.minLeafSize_ << token::SPACE
+            << tree.boundBoxes_ << token::SPACE
+            << tree.addressing_ << token::SPACE;
+    }
+    else
+    {
+        os.write
+        (
+            reinterpret_cast<const char*>(&tree.maxLevel_),
+            sizeof(tree.maxLevel_)
+          + sizeof(tree.minLeafSize_)
+        );
+        os  << tree.boundBoxes_
+            << tree.addressing_;
+    }
+
+    os.check(FUNCTION_NAME);
+    return os;
+}
+
+
+template<class Type>
+Foam::Istream& Foam::operator>>(Istream& is, AABBTree<Type>& tree)
+{
+    if (is.format() == IOstreamOption::ASCII)
+    {
+        is  >> tree.maxLevel_
+            >> tree.minLeafSize_;
+    }
+    else
+    {
+        is.beginRawRead();
+
+        readRawLabel(is, &tree.maxLevel_);
+        readRawLabel(is, &tree.minLeafSize_);
+
+        is.endRawRead();
+    }
+
+    is  >> tree.boundBoxes_
+        >> tree.addressing_;
+
+    is.check(FUNCTION_NAME);
+    return is;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/algorithms/AABBTree/AABBTree.H

@@ -0,0 +1,198 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::AABBTree
+
+Description
+    Templated tree of axis-aligned bounding boxes (AABB)
+
+    Designed to be templated on either faces or cells, the AABBTree will
+    decompose the input into a tree of AABB's.  The maximum number of tree
+    levels and minimum number of objects per leaf are provided on construction,
+    and the contents (addressing) is stored.
+
+SourceFiles
+    AABBTree.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef AABBTree_H
+#define AABBTree_H
+
+#include "labelList.H"
+#include "labelPair.H"
+#include "DynamicList.H"
+#include "pointField.H"
+#include "treeBoundBox.H"
+#include "Ostream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class Type>
+class AABBTree;
+
+template<class Type>
+Istream& operator>>(Istream&, AABBTree<Type>&);
+
+template<class Type>
+Ostream& operator<<(Ostream&, const AABBTree<Type>&);
+
+/*---------------------------------------------------------------------------*\
+                          Class AABBTree Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class AABBTree
+{
+
+protected:
+
+    // Protected Data
+
+        //- Tolerance
+        static scalar tolerance_;
+
+        //- Maximum tree level
+        label maxLevel_;
+
+        //- Minimum points per leaf
+        label minLeafSize_;
+
+        //- Bounding boxes making up the tree
+        List<treeBoundBox> boundBoxes_;
+
+        //- Leaf addressing
+        List<labelList> addressing_;
+
+
+    // Protected Member Functions
+
+        //- Write OBJ file of bounding box
+        void writeOBJ
+        (
+            const bool writeLinesOnly,
+            const treeBoundBox& bb,
+            label& vertI,
+            Ostream& os
+        ) const;
+
+        //- Write OBJ for all bounding boxes
+        void writeOBJ
+        (
+            const bool leavesOnly,
+            const bool writeLinesOnly,
+            const treeBoundBox& bb,
+            const label nodeI,
+            const List<Pair<treeBoundBox>>& bbs,
+            const List<Pair<label>>& nodes,
+            label& vertI,
+            Ostream& os
+        ) const;
+
+        //- Create the bounding boxes by interrogating points
+        void createBoxes
+        (
+            const bool equalBinSize,
+            const label level,
+            const List<Type>& objects,
+            const pointField& points,
+            const DynamicList<label>& objectIDs,
+            const treeBoundBox& bb,
+            const label nodeI,
+
+            DynamicList<Pair<treeBoundBox>>& bbs,
+            DynamicList<labelPair>& nodes,
+            DynamicList<labelList>& addressing
+        ) const;
+
+
+public:
+
+    // Constructors
+
+        //- Null constructor
+        AABBTree();
+
+        //- Construct from components
+        //  equalBinSize: divide into equal number of elements or equal span
+        AABBTree
+        (
+            const UList<Type>& objects,
+            const pointField& points,
+            const bool equalBinSize = true,
+            const label maxLevel = 3,
+            const label minBinSize = 100
+        );
+
+
+    // Public Member Functions
+
+        // Access
+
+            //- Return the bounding boxes making up the tree
+            const List<treeBoundBox>& boundBoxes() const;
+
+            //- Return the contents addressing
+            const List<labelList>& addressing() const;
+
+
+        // Evaluation
+
+            //- Determine whether a point is inside the bounding boxes
+            bool pointInside(const point& pt) const;
+
+            //- Determine whether a bounding box overlaps the tree bounding
+            //- boxes
+            bool overlaps(const boundBox& bbIn) const;
+
+
+    // IOstream operators
+
+        friend Istream& operator>> <Type>(Istream&, AABBTree&);
+        friend Ostream& operator<< <Type>(Ostream&, const AABBTree&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "AABBTree.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.C

@@ -0,0 +1,273 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "treeDataEdge.H"
+#include "indexedOctree.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(treeDataEdge, 0);
+}
+
+
+// * * * * * * * * * * * * * 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)
+        {
+            bbs_[i] = calcBb(edgeLabels_[i]);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::treeDataEdge::treeDataEdge
+(
+    const bool cacheBb,
+    const edgeList& edges,
+    const pointField& points,
+    const labelUList& edgeLabels
+)
+:
+    edges_(edges),
+    points_(points),
+    edgeLabels_(edgeLabels),
+    cacheBb_(cacheBb)
+{
+    update();
+}
+
+
+Foam::treeDataEdge::treeDataEdge
+(
+    const bool cacheBb,
+    const edgeList& edges,
+    const pointField& points,
+    labelList&& edgeLabels
+)
+:
+    edges_(edges),
+    points_(points),
+    edgeLabels_(std::move(edgeLabels)),
+    cacheBb_(cacheBb)
+{
+    update();
+}
+
+
+Foam::treeDataEdge::findNearestOp::findNearestOp
+(
+    const indexedOctree<treeDataEdge>& tree
+)
+:
+    tree_(tree)
+{}
+
+
+Foam::treeDataEdge::findIntersectOp::findIntersectOp
+(
+    const indexedOctree<treeDataEdge>& tree
+)
+{}
+
+
+// * * * * * * * * * * * * * * * 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<treeDataEdge>& 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()
+(
+    const labelUList& indices,
+    const point& sample,
+
+    scalar& nearestDistSqr,
+    label& minIndex,
+    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);
+
+        const scalar distSqr = sqr(nearHit.distance());
+
+        if (distSqr < nearestDistSqr)
+        {
+            nearestDistSqr = distSqr;
+            minIndex = index;
+            nearestPoint = nearHit.rawPoint();
+        }
+    }
+}
+
+
+void Foam::treeDataEdge::findNearestOp::operator()
+(
+    const labelUList& indices,
+    const linePointRef& ln,
+
+    treeBoundBox& tightest,
+    label& minIndex,
+    point& linePoint,
+    point& nearestPoint
+) const
+{
+    const treeDataEdge& shape = tree_.shapes();
+
+    // Best so far
+    scalar nearestDistSqr = magSqr(linePoint - 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);
+
+        if (distSqr < nearestDistSqr)
+        {
+            nearestDistSqr = distSqr;
+            minIndex = index;
+            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;
+            }
+        }
+    }
+}
+
+
+bool Foam::treeDataEdge::findIntersectOp::operator()
+(
+    const label index,
+    const point& start,
+    const point& end,
+    point& result
+) const
+{
+    NotImplemented;
+    return false;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.H

@@ -0,0 +1,235 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::treeDataEdge
+
+Description
+    Holds data for octree to work on an edges subset.
+
+SourceFiles
+    treeDataEdge.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_treeDataEdge_H
+#define Foam_treeDataEdge_H
+
+#include "treeBoundBoxList.H"
+#include "line.H"
+#include "volumeType.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+template<class Type> class indexedOctree;
+
+/*---------------------------------------------------------------------------*\
+                           Class treeDataEdge Declaration
+\*---------------------------------------------------------------------------*/
+
+class treeDataEdge
+{
+    // Static Data
+
+        //- Tolerance on linear dimensions
+        static scalar tol;
+
+
+    // Private Data
+
+        //- Reference to edgeList
+        const edgeList& edges_;
+
+        //- Reference to points
+        const pointField& points_;
+
+        //- Labels of edges
+        const labelList edgeLabels_;
+
+        //- Whether to precalculate and store face bounding box
+        const bool cacheBb_;
+
+        //- Bbs for all above edges (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:
+
+
+    class findNearestOp
+    {
+        const indexedOctree<treeDataEdge>& tree_;
+
+    public:
+
+        findNearestOp(const indexedOctree<treeDataEdge>& tree);
+
+        void operator()
+        (
+            const labelUList& indices,
+            const point& sample,
+
+            scalar& nearestDistSqr,
+            label& minIndex,
+            point& nearestPoint
+        ) const;
+
+        void operator()
+        (
+            const labelUList& indices,
+            const linePointRef& ln,
+
+            treeBoundBox& tightest,
+            label& minIndex,
+            point& linePoint,
+            point& nearestPoint
+        ) const;
+    };
+
+
+    class findIntersectOp
+    {
+    public:
+
+        findIntersectOp(const indexedOctree<treeDataEdge>& tree);
+
+        //- Calculate intersection of triangle with ray. Sets result
+        //  accordingly
+        bool operator()
+        (
+            const label index,
+            const point& start,
+            const point& end,
+            point& intersectionPoint
+        ) const;
+    };
+
+
+    // Declare name of the class and its debug switch
+    ClassName("treeDataEdge");
+
+
+    // Constructors
+
+        //- Construct from selected edges.
+        //  \note Holds references to edges and points!
+        treeDataEdge
+        (
+            const bool cacheBb,
+            const edgeList& edges,
+            const pointField& points,
+            const labelUList& edgeLabels
+        );
+
+        //- Construct from selected edges, transferring contents.
+        //  \note Holds references to edges and points!
+        treeDataEdge
+        (
+            const bool cacheBb,
+            const edgeList& edges,
+            const pointField& points,
+            labelList&& edgeLabels
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            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;
+
+
+        // Search
+
+            //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+            //  Only makes sense for closed surfaces.
+            volumeType getVolumeType
+            (
+                const indexedOctree<treeDataEdge>&,
+                const point&
+            ) const;
+
+            //- Does (bb of) shape at index overlap bb
+            bool overlaps
+            (
+                const label index,
+                const treeBoundBox& sampleBb
+            ) const;
+
+            //- Does (bb of) shape at index overlap bb
+            bool overlaps
+            (
+                const label index,
+                const point& centre,
+                const scalar radiusSqr
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.C

@@ -0,0 +1,233 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "treeDataPoint.H"
+#include "indexedOctree.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(treeDataPoint, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::treeDataPoint::treeDataPoint(const pointField& points)
+:
+    points_(points),
+    useSubset_(false)
+{}
+
+
+Foam::treeDataPoint::treeDataPoint
+(
+    const pointField& points,
+    const labelUList& pointLabels,
+    const bool useSubsetPoints
+)
+:
+    points_(points),
+    pointLabels_(pointLabels),
+    useSubset_(useSubsetPoints)
+{}
+
+
+Foam::treeDataPoint::treeDataPoint
+(
+    const pointField& points,
+    labelList&& pointLabels,
+    const bool useSubsetPoints
+)
+:
+    points_(points),
+    pointLabels_(std::move(pointLabels)),
+    useSubset_(useSubsetPoints)
+{}
+
+
+Foam::treeDataPoint::findNearestOp::findNearestOp
+(
+    const indexedOctree<treeDataPoint>& tree
+)
+:
+    tree_(tree)
+{}
+
+
+Foam::treeDataPoint::findIntersectOp::findIntersectOp
+(
+    const indexedOctree<treeDataPoint>& tree
+)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::pointField Foam::treeDataPoint::shapePoints() const
+{
+    if (useSubset_)
+    {
+        return pointField(points_, pointLabels_);
+    }
+
+    return points_;
+}
+
+
+Foam::volumeType Foam::treeDataPoint::getVolumeType
+(
+    const indexedOctree<treeDataPoint>& oc,
+    const point& sample
+) const
+{
+    return volumeType::UNKNOWN;
+}
+
+
+bool Foam::treeDataPoint::overlaps
+(
+    const label index,
+    const treeBoundBox& cubeBb
+) const
+{
+    return cubeBb.contains(shapePoint(index));
+}
+
+
+bool Foam::treeDataPoint::overlaps
+(
+    const label index,
+    const point& centre,
+    const scalar radiusSqr
+) const
+{
+    return (magSqr(shapePoint(index) - centre) <= radiusSqr);
+}
+
+
+void Foam::treeDataPoint::findNearestOp::operator()
+(
+    const labelUList& indices,
+    const point& sample,
+
+    scalar& nearestDistSqr,
+    label& minIndex,
+    point& nearestPoint
+) const
+{
+    const treeDataPoint& shape = tree_.shapes();
+
+    for (const label index : indices)
+    {
+        const point& pt = shape.shapePoint(index);
+
+        const scalar distSqr = magSqr(pt - sample);
+
+        if (distSqr < nearestDistSqr)
+        {
+            nearestDistSqr = distSqr;
+            minIndex = index;
+            nearestPoint = pt;
+        }
+    }
+}
+
+
+void Foam::treeDataPoint::findNearestOp::operator()
+(
+    const labelUList& indices,
+    const linePointRef& ln,
+
+    treeBoundBox& tightest,
+    label& minIndex,
+    point& linePoint,
+    point& nearestPoint
+) const
+{
+    const treeDataPoint& shape = tree_.shapes();
+
+    // Best so far
+    scalar nearestDistSqr = GREAT;
+    if (minIndex >= 0)
+    {
+        nearestDistSqr = magSqr(linePoint - nearestPoint);
+    }
+
+    for (const label index : indices)
+    {
+        const point& shapePt = shape.shapePoint(index);
+
+        if (tightest.contains(shapePt))
+        {
+            // Nearest point on line
+            pointHit pHit = ln.nearestDist(shapePt);
+            const scalar distSqr = sqr(pHit.distance());
+
+            if (distSqr < nearestDistSqr)
+            {
+                nearestDistSqr = distSqr;
+                minIndex = index;
+                linePoint = pHit.rawPoint();
+                nearestPoint = shapePt;
+
+                {
+                    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();
+                }
+            }
+        }
+    }
+}
+
+
+bool Foam::treeDataPoint::findIntersectOp::operator()
+(
+    const label index,
+    const point& start,
+    const point& end,
+    point& result
+) const
+{
+    NotImplemented;
+    return false;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.H

@@ -0,0 +1,271 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2013 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::treeDataPoint
+
+Description
+    Holds (reference to) pointField. Encapsulation of data needed for
+    octree searches.
+    Used for searching for nearest point. No bounding boxes around points.
+    Only overlaps and calcNearest are implemented, rest makes little sense.
+
+    Optionally works on subset of points.
+
+SourceFiles
+    treeDataPoint.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_treeDataPoint_H
+#define Foam_treeDataPoint_H
+
+#include "pointField.H"
+#include "treeBoundBox.H"
+#include "line.H"
+#include "volumeType.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+template<class Type> class indexedOctree;
+
+/*---------------------------------------------------------------------------*\
+                        Class treeDataPoint Declaration
+\*---------------------------------------------------------------------------*/
+
+class treeDataPoint
+{
+    // Private Data
+
+        const pointField& points_;
+
+        //- Subset of points to work on (or empty)
+        const labelList pointLabels_;
+
+        const bool useSubset_;
+
+public:
+
+
+    class findNearestOp
+    {
+        const indexedOctree<treeDataPoint>& tree_;
+
+    public:
+
+        findNearestOp(const indexedOctree<treeDataPoint>& tree);
+
+        void operator()
+        (
+            const labelUList& indices,
+            const point& sample,
+
+            scalar& nearestDistSqr,
+            label& minIndex,
+            point& nearestPoint
+        ) const;
+
+        void operator()
+        (
+            const labelUList& indices,
+            const linePointRef& ln,
+
+            treeBoundBox& tightest,
+            label& minIndex,
+            point& linePoint,
+            point& nearestPoint
+        ) const;
+    };
+
+
+    class findIntersectOp
+    {
+    public:
+
+        findIntersectOp(const indexedOctree<treeDataPoint>& tree);
+
+        //- Calculate intersection of triangle with ray.
+        //  Sets result accordingly
+        bool operator()
+        (
+            const label index,
+            const point& start,
+            const point& end,
+            point& intersectionPoint
+        ) const;
+    };
+
+
+    // Declare name of the class and its debug switch
+    ClassName("treeDataPoint");
+
+
+    // Constructors
+
+        //- Construct from pointField
+        //  \note Holds reference to the points!
+        explicit treeDataPoint(const pointField& points);
+
+        //- Construct from subset of pointField, copies point ids
+        //  \note Holds reference to the points!
+        treeDataPoint
+        (
+            const pointField& points,
+            const labelUList& pointLabels,
+            const bool useSubsetPoints = true
+        );
+
+        //- Construct from subset of pointField, moves point ids
+        //  \note Holds reference to the points!
+        treeDataPoint
+        (
+            const pointField& points,
+            labelList&& pointLabels,
+            const bool useSubsetPoints = true
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            //- 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;
+
+
+        // Search
+
+            //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+            //  Only makes sense for closed surfaces.
+            volumeType getVolumeType
+            (
+                const indexedOctree<treeDataPoint>& os,
+                const point& sample
+            ) const;
+
+            //- Does (bb of) shape at index overlap bb
+            bool overlaps
+            (
+                const label index,
+                const treeBoundBox& sampleBb
+            ) const;
+
+            //- Does shape at index overlap the sphere
+            bool overlaps
+            (
+                const label index,
+                const point& centre,
+                const scalar radiusSqr
+            ) const;
+
+
+    // Member Operators
+
+        //- The point at the specified index
+        inline const point& operator[](const label index) const
+        {
+            return shapePoint(index);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //