src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "treeDataCell.H"
  27 #include "indexedOctree.H"
  28 #include "primitiveMesh.H"
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 defineTypeNameAndDebug(Foam::treeDataCell, 0);
  33
  34
  35 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  36
  37 Foam::treeBoundBox Foam::treeDataCell::calcCellBb(const label cellI) const
  38 {
  39     const cellList& cells = mesh_.cells();
  40     const faceList& faces = mesh_.faces();
  41     const pointField& points = mesh_.points();
  42
  43     treeBoundBox cellBb
  44     (
  45         vector(GREAT, GREAT, GREAT),
  46         vector(-GREAT, -GREAT, -GREAT)
  47     );
  48
  49     const cell& cFaces = cells[cellI];
  50
  51     forAll(cFaces, cFaceI)
  52     {
  53         const face& f = faces[cFaces[cFaceI]];
  54
  55         forAll(f, fp)
  56         {
  57             const point& p = points[f[fp]];
  58
  59             cellBb.min() = min(cellBb.min(), p);
  60             cellBb.max() = max(cellBb.max(), p);
  61         }
  62     }
  63     return cellBb;
  64 }
  65
  66
  67 void Foam::treeDataCell::update()
  68 {
  69     if (cacheBb_)
  70     {
  71         bbs_.setSize(cellLabels_.size());
  72
  73         forAll(cellLabels_, i)
  74         {
  75             bbs_[i] = calcCellBb(cellLabels_[i]);
  76         }
  77     }
  78 }
  79
  80
  81 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  82
  83 Foam::treeDataCell::treeDataCell
  84 (
  85     const bool cacheBb,
  86     const primitiveMesh& mesh,
  87     const labelUList& cellLabels
  88 )
  89 :
  90     mesh_(mesh),
  91     cellLabels_(cellLabels),
  92     cacheBb_(cacheBb)
  93 {
  94     update();
  95 }
  96
  97
  98 Foam::treeDataCell::treeDataCell
  99 (
 100     const bool cacheBb,
 101     const primitiveMesh& mesh,
 102     const Xfer<labelList>& cellLabels
 103 )
 104 :
 105     mesh_(mesh),
 106     cellLabels_(cellLabels),
 107     cacheBb_(cacheBb)
 108 {
 109     update();
 110 }
 111
 112
 113 Foam::treeDataCell::treeDataCell
 114 (
 115     const bool cacheBb,
 116     const primitiveMesh& mesh
 117 )
 118 :
 119     mesh_(mesh),
 120     cellLabels_(identity(mesh_.nCells())),
 121     cacheBb_(cacheBb)
 122 {
 123     update();
 124 }
 125
 126
 127 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 128
 129 Foam::pointField Foam::treeDataCell::shapePoints() const
 130 {
 131     pointField cc(cellLabels_.size());
 132
 133     forAll(cellLabels_, i)
 134     {
 135         cc[i] = mesh_.cellCentres()[cellLabels_[i]];
 136     }
 137
 138     return cc;
 139 }
 140
 141
 142 bool Foam::treeDataCell::overlaps
 143 (
 144     const label index,
 145     const treeBoundBox& cubeBb
 146 ) const
 147 {
 148     if (cacheBb_)
 149     {
 150         return cubeBb.overlaps(bbs_[index]);
 151     }
 152     else
 153     {
 154         return cubeBb.overlaps(calcCellBb(cellLabels_[index]));
 155     }
 156 }
 157
 158
 159 bool Foam::treeDataCell::contains
 160 (
 161     const label index,
 162     const point& sample
 163 ) const
 164 {
 165     return mesh_.pointInCell(sample, cellLabels_[index]);
 166 }
 167
 168
 169 void Foam::treeDataCell::findNearest
 170 (
 171     const labelUList& indices,
 172     const point& sample,
 173
 174     scalar& nearestDistSqr,
 175     label& minIndex,
 176     point& nearestPoint
 177 ) const
 178 {
 179     forAll(indices, i)
 180     {
 181         label index = indices[i];
 182         label cellI = cellLabels_[index];
 183         scalar distSqr = magSqr(sample - mesh_.cellCentres()[cellI]);
 184
 185         if (distSqr < nearestDistSqr)
 186         {
 187             nearestDistSqr = distSqr;
 188             minIndex = index;
 189             nearestPoint = mesh_.cellCentres()[cellI];
 190         }
 191     }
 192 }
 193
 194
 195 bool Foam::treeDataCell::intersects
 196 (
 197     const label index,
 198     const point& start,
 199     const point& end,
 200     point& intersectionPoint
 201 ) const
 202 {
 203     // Do quick rejection test
 204     if (cacheBb_)
 205     {
 206         const treeBoundBox& cellBb = bbs_[index];
 207
 208         if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0)
 209         {
 210             // start and end in same block outside of cellBb.
 211             return false;
 212         }
 213     }
 214     else
 215     {
 216         const treeBoundBox cellBb = calcCellBb(cellLabels_[index]);
 217
 218         if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0)
 219         {
 220             // start and end in same block outside of cellBb.
 221             return false;
 222         }
 223     }
 224
 225
 226     // Do intersection with all faces of cell
 227     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 228
 229     // Disable picking up intersections behind us.
 230     scalar oldTol = intersection::setPlanarTol(0.0);
 231
 232     const cell& cFaces = mesh_.cells()[cellLabels_[index]];
 233
 234     const vector dir(end - start);
 235     scalar minDistSqr = magSqr(dir);
 236     bool hasMin = false;
 237
 238     forAll(cFaces, i)
 239     {
 240         const face& f = mesh_.faces()[cFaces[i]];
 241
 242         pointHit inter = f.ray
 243         (
 244             start,
 245             dir,
 246             mesh_.points(),
 247             intersection::HALF_RAY
 248         );
 249
 250         if (inter.hit() && sqr(inter.distance()) <= minDistSqr)
 251         {
 252             // Note: no extra test on whether intersection is in front of us
 253             // since using half_ray AND zero tolerance. (note that tolerance
 254             // is used to look behind us)
 255             minDistSqr = sqr(inter.distance());
 256             intersectionPoint = inter.hitPoint();
 257             hasMin = true;
 258         }
 259     }
 260
 261     // Restore picking tolerance
 262     intersection::setPlanarTol(oldTol);
 263
 264     return hasMin;
 265 }
 266
 267
 268 // ************************************************************************* //