/*---------------------------------------------------------------------------*\ ========= | \\ / 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) 2017-2020 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 "treeBoundBox.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // const Foam::faceList Foam::treeBoundBox::faces ({ face({0, 4, 6, 2}), // 0: x-min, left face({1, 3, 7, 5}), // 1: x-max, right face({0, 1, 5, 4}), // 2: y-min, bottom face({2, 6, 7, 3}), // 3: y-max, top face({0, 2, 3, 1}), // 4: z-min, back face({4, 5, 7, 6}) // 5: z-max, front }); const Foam::edgeList Foam::treeBoundBox::edges ({ {0, 1}, // 0 {1, 3}, {2, 3}, // 2 {0, 2}, {4, 5}, // 4 {5, 7}, {6, 7}, // 6 {4, 6}, {0, 4}, // 8 {1, 5}, {3, 7}, // 10 {2, 6} }); // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::treeBoundBox::treeBoundBox(const UList<point>& points) : boundBox(points, false) { if (points.empty()) { WarningInFunction << "No bounding box for zero-sized pointField" << nl; } } Foam::treeBoundBox::treeBoundBox ( const UList<point>& points, const labelUList& indices ) : boundBox(points, indices, false) { if (points.empty() || indices.empty()) { WarningInFunction << "No bounding box for zero-sized pointField" << nl; } } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // Foam::tmp<Foam::pointField> Foam::treeBoundBox::points() const { auto tpts = tmp<pointField>::New(8); auto& pts = tpts.ref(); forAll(pts, octant) { pts[octant] = corner(octant); } return tpts; } Foam::treeBoundBox Foam::treeBoundBox::subBbox(const direction octant) const { return subBbox(centre(), octant); } Foam::treeBoundBox Foam::treeBoundBox::subBbox ( const point& mid, const direction octant ) const { if (octant > 7) { FatalErrorInFunction << "octant should be [0..7]" << abort(FatalError); } // start with a copy of this bounding box and adjust limits accordingly treeBoundBox subBb(*this); point& bbMin = subBb.min(); point& bbMax = subBb.max(); if (octant & treeBoundBox::RIGHTHALF) { bbMin.x() = mid.x(); // mid -> max } else { bbMax.x() = mid.x(); // min -> mid } if (octant & treeBoundBox::TOPHALF) { bbMin.y() = mid.y(); // mid -> max } else { bbMax.y() = mid.y(); // min -> mid } if (octant & treeBoundBox::FRONTHALF) { bbMin.z() = mid.z(); // mid -> max } else { bbMax.z() = mid.z(); // min -> mid } return subBb; } bool Foam::treeBoundBox::intersects ( const point& overallStart, const vector& overallVec, const point& start, const point& end, point& pt, direction& ptOnFaces ) const { // Sutherlands algorithm: // loop // - start = intersection of line with one of the planes bounding // the bounding box // - stop if start inside bb (return true) // - stop if start and end in same 'half' (e.g. both above bb) // (return false) // // Uses posBits to efficiently determine 'half' in which start and end // point are. // // Note: // - sets coordinate to exact position: e.g. pt.x() = min().x() // since plane intersect routine might have truncation error. // This makes sure that posBits tests 'inside' const direction endBits = posBits(end); pt = start; // Allow maximum of 3 clips. for (label i = 0; i < 4; ++i) { direction ptBits = posBits(pt); if (ptBits == 0) { // pt inside bb ptOnFaces = faceBits(pt); return true; } if ((ptBits & endBits) != 0) { // pt and end in same block outside of bb ptOnFaces = faceBits(pt); return false; } if (ptBits & LEFTBIT) { // Intersect with plane V=min, n=-1,0,0 if (Foam::mag(overallVec.x()) > VSMALL) { scalar s = (min().x() - overallStart.x())/overallVec.x(); pt.x() = min().x(); pt.y() = overallStart.y() + overallVec.y()*s; pt.z() = overallStart.z() + overallVec.z()*s; } else { // Vector not in x-direction. But still intersecting bb planes. // So must be close - just snap to bb. pt.x() = min().x(); } } else if (ptBits & RIGHTBIT) { // Intersect with plane V=max, n=1,0,0 if (Foam::mag(overallVec.x()) > VSMALL) { scalar s = (max().x() - overallStart.x())/overallVec.x(); pt.x() = max().x(); pt.y() = overallStart.y() + overallVec.y()*s; pt.z() = overallStart.z() + overallVec.z()*s; } else { pt.x() = max().x(); } } else if (ptBits & BOTTOMBIT) { // Intersect with plane V=min, n=0,-1,0 if (Foam::mag(overallVec.y()) > VSMALL) { scalar s = (min().y() - overallStart.y())/overallVec.y(); pt.x() = overallStart.x() + overallVec.x()*s; pt.y() = min().y(); pt.z() = overallStart.z() + overallVec.z()*s; } else { pt.x() = min().y(); } } else if (ptBits & TOPBIT) { // Intersect with plane V=max, n=0,1,0 if (Foam::mag(overallVec.y()) > VSMALL) { scalar s = (max().y() - overallStart.y())/overallVec.y(); pt.x() = overallStart.x() + overallVec.x()*s; pt.y() = max().y(); pt.z() = overallStart.z() + overallVec.z()*s; } else { pt.y() = max().y(); } } else if (ptBits & BACKBIT) { // Intersect with plane V=min, n=0,0,-1 if (Foam::mag(overallVec.z()) > VSMALL) { scalar s = (min().z() - overallStart.z())/overallVec.z(); pt.x() = overallStart.x() + overallVec.x()*s; pt.y() = overallStart.y() + overallVec.y()*s; pt.z() = min().z(); } else { pt.z() = min().z(); } } else if (ptBits & FRONTBIT) { // Intersect with plane V=max, n=0,0,1 if (Foam::mag(overallVec.z()) > VSMALL) { scalar s = (max().z() - overallStart.z())/overallVec.z(); pt.x() = overallStart.x() + overallVec.x()*s; pt.y() = overallStart.y() + overallVec.y()*s; pt.z() = max().z(); } else { pt.z() = max().z(); } } } // Can end up here if the end point is on the edge of the boundBox return true; } bool Foam::treeBoundBox::intersects ( const point& start, const point& end, point& pt ) const { direction ptBits; return intersects(start, end-start, start, end, pt, ptBits); } bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const { // Compare all components against min and max of bb for (direction cmpt=0; cmpt < point::nComponents; ++cmpt) { if (pt[cmpt] < min()[cmpt]) { return false; } else if (pt[cmpt] == min()[cmpt]) { // On edge. Outside if direction points outwards. if (dir[cmpt] < 0) { return false; } } if (pt[cmpt] > max()[cmpt]) { return false; } else if (pt[cmpt] == max()[cmpt]) { // On edge. Outside if direction points outwards. if (dir[cmpt] > 0) { return false; } } } // All components inside bb return true; } Foam::direction Foam::treeBoundBox::faceBits(const point& pt) const { direction octant = 0; if (pt.x() == min().x()) { octant |= LEFTBIT; } else if (pt.x() == max().x()) { octant |= RIGHTBIT; } if (pt.y() == min().y()) { octant |= BOTTOMBIT; } else if (pt.y() == max().y()) { octant |= TOPBIT; } if (pt.z() == min().z()) { octant |= BACKBIT; } else if (pt.z() == max().z()) { octant |= FRONTBIT; } return octant; } Foam::direction Foam::treeBoundBox::posBits(const point& pt) const { direction octant = 0; if (pt.x() < min().x()) { octant |= LEFTBIT; } else if (pt.x() > max().x()) { octant |= RIGHTBIT; } if (pt.y() < min().y()) { octant |= BOTTOMBIT; } else if (pt.y() > max().y()) { octant |= TOPBIT; } if (pt.z() < min().z()) { octant |= BACKBIT; } else if (pt.z() > max().z()) { octant |= FRONTBIT; } return octant; } void Foam::treeBoundBox::calcExtremities ( const point& pt, point& nearest, point& furthest ) const { for (direction cmpt=0; cmpt < point::nComponents; ++cmpt) { if ( Foam::mag(min()[cmpt] - pt[cmpt]) < Foam::mag(max()[cmpt] - pt[cmpt]) ) { nearest[cmpt] = min()[cmpt]; furthest[cmpt] = max()[cmpt]; } else { nearest[cmpt] = max()[cmpt]; furthest[cmpt] = min()[cmpt]; } } } Foam::scalar Foam::treeBoundBox::maxDist(const point& pt) const { point near, far; calcExtremities(pt, near, far); return Foam::mag(far - pt); } Foam::label Foam::treeBoundBox::distanceCmp ( const point& pt, const treeBoundBox& other ) const { // // Distance point <-> nearest and furthest away vertex of this // point nearThis, farThis; // get nearest and furthest away vertex calcExtremities(pt, nearThis, farThis); const scalar minDistThis = sqr(nearThis.x() - pt.x()) + sqr(nearThis.y() - pt.y()) + sqr(nearThis.z() - pt.z()); const scalar maxDistThis = sqr(farThis.x() - pt.x()) + sqr(farThis.y() - pt.y()) + sqr(farThis.z() - pt.z()); // // Distance point <-> other // point nearOther, farOther; // get nearest and furthest away vertex other.calcExtremities(pt, nearOther, farOther); const scalar minDistOther = sqr(nearOther.x() - pt.x()) + sqr(nearOther.y() - pt.y()) + sqr(nearOther.z() - pt.z()); const scalar maxDistOther = sqr(farOther.x() - pt.x()) + sqr(farOther.y() - pt.y()) + sqr(farOther.z() - pt.z()); // // Categorize // if (maxDistThis < minDistOther) { // All vertices of this are nearer to point than any vertex of other return -1; } else if (minDistThis > maxDistOther) { // All vertices of this are further from point than any vertex of other return 1; } else { // Mixed bag return 0; } } // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // Foam::Ostream& Foam::operator<<(Ostream& os, const treeBoundBox& bb) { return os << static_cast<const boundBox&>(bb); } Foam::Istream& Foam::operator>>(Istream& is, treeBoundBox& bb) { return is >> static_cast<boundBox&>(bb); } // ************************************************************************* //