src/meshTools/indexedOctree/treeDataCell.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "treeDataCell.H"
  28 #include "indexedOctree.H"
  29 #include <OpenFOAM/primitiveMesh.H>
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 defineTypeNameAndDebug(Foam::treeDataCell, 0);
  34
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38 Foam::treeBoundBox Foam::treeDataCell::calcCellBb(const label cellI) const
  39 {
  40     const cellList& cells = mesh_.cells();
  41     const faceList& faces = mesh_.faces();
  42     const pointField& points = mesh_.points();
  43
  44     treeBoundBox cellBb
  45     (
  46         vector(GREAT, GREAT, GREAT),
  47         vector(-GREAT, -GREAT, -GREAT)
  48     );
  49
  50     const cell& cFaces = cells[cellI];
  51
  52     forAll(cFaces, cFaceI)
  53     {
  54         const face& f = faces[cFaces[cFaceI]];
  55
  56         forAll(f, fp)
  57         {
  58             const point& p = points[f[fp]];
  59
  60             cellBb.min() = min(cellBb.min(), p);
  61             cellBb.max() = max(cellBb.max(), p);
  62         }
  63     }
  64     return cellBb;
  65 }
  66
  67
  68 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  69
  70 // Construct from components
  71 Foam::treeDataCell::treeDataCell
  72 (
  73     const bool cacheBb,
  74     const primitiveMesh& mesh,
  75     const labelList& cellLabels
  76 )
  77 :
  78     mesh_(mesh),
  79     cellLabels_(cellLabels),
  80     cacheBb_(cacheBb)
  81 {
  82     if (cacheBb_)
  83     {
  84         bbs_.setSize(cellLabels_.size());
  85
  86         forAll(cellLabels_, i)
  87         {
  88             bbs_[i] = calcCellBb(cellLabels_[i]);
  89         }
  90     }
  91 }
  92
  93
  94 Foam::treeDataCell::treeDataCell
  95 (
  96     const bool cacheBb,
  97     const primitiveMesh& mesh
  98 )
  99 :
 100     mesh_(mesh),
 101     cellLabels_(identity(mesh_.nCells())),
 102     cacheBb_(cacheBb)
 103 {
 104     if (cacheBb_)
 105     {
 106         bbs_.setSize(cellLabels_.size());
 107
 108         forAll(cellLabels_, i)
 109         {
 110             bbs_[i] = calcCellBb(cellLabels_[i]);
 111         }
 112     }
 113 }
 114
 115
 116 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 117
 118 Foam::pointField Foam::treeDataCell::points() const
 119 {
 120     pointField cc(cellLabels_.size());
 121
 122     forAll(cellLabels_, i)
 123     {
 124         cc[i] = mesh_.cellCentres()[cellLabels_[i]];
 125     }
 126
 127     return cc;
 128 }
 129
 130
 131 // Check if any point on shape is inside cubeBb.
 132 bool Foam::treeDataCell::overlaps
 133 (
 134     const label index,
 135     const treeBoundBox& cubeBb
 136 ) const
 137 {
 138     if (cacheBb_)
 139     {
 140         return cubeBb.overlaps(bbs_[index]);
 141     }
 142     else
 143     {
 144         return cubeBb.overlaps(calcCellBb(cellLabels_[index]));
 145     }
 146 }
 147
 148
 149 // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
 150 // nearestPoint.
 151 void Foam::treeDataCell::findNearest
 152 (
 153     const labelList& indices,
 154     const point& sample,
 155
 156     scalar& nearestDistSqr,
 157     label& minIndex,
 158     point& nearestPoint
 159 ) const
 160 {
 161     forAll(indices, i)
 162     {
 163         label index = indices[i];
 164         label cellI = cellLabels_[index];
 165         scalar distSqr = magSqr(sample - mesh_.cellCentres()[cellI]);
 166
 167         if (distSqr < nearestDistSqr)
 168         {
 169             nearestDistSqr = distSqr;
 170             minIndex = index;
 171             nearestPoint = mesh_.cellCentres()[cellI];
 172         }
 173     }
 174 }
 175
 176
 177 bool Foam::treeDataCell::intersects
 178 (
 179     const label index,
 180     const point& start,
 181     const point& end,
 182     point& intersectionPoint
 183 ) const
 184 {
 185     // Do quick rejection test
 186     if (cacheBb_)
 187     {
 188         const treeBoundBox& cellBb = bbs_[index];
 189
 190         if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0)
 191         {
 192             // start and end in same block outside of cellBb.
 193             return false;
 194         }
 195     }
 196     else
 197     {
 198         const treeBoundBox cellBb = calcCellBb(cellLabels_[index]);
 199
 200         if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0)
 201         {
 202             // start and end in same block outside of cellBb.
 203             return false;
 204         }
 205     }
 206
 207
 208     // Do intersection with all faces of cell
 209     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 210
 211     // Disable picking up intersections behind us.
 212     scalar oldTol = intersection::setPlanarTol(0.0);
 213
 214     const cell& cFaces = mesh_.cells()[cellLabels_[index]];
 215
 216     const vector dir(end - start);
 217     scalar minDistSqr = magSqr(dir);
 218     bool hasMin = false;
 219
 220     forAll(cFaces, i)
 221     {
 222         const face& f = mesh_.faces()[cFaces[i]];
 223
 224         pointHit inter = f.ray
 225         (
 226             start,
 227             dir,
 228             mesh_.points(),
 229             intersection::HALF_RAY
 230         );
 231
 232         if (inter.hit() && sqr(inter.distance()) <= minDistSqr)
 233         {
 234             // Note: no extra test on whether intersection is in front of us
 235             // since using half_ray AND zero tolerance. (note that tolerance
 236             // is used to look behind us)
 237             minDistSqr = sqr(inter.distance());
 238             intersectionPoint = inter.hitPoint();
 239             hasMin = true;
 240         }
 241     }
 242
 243     // Restore picking tolerance
 244     intersection::setPlanarTol(oldTol);
 245
 246     return hasMin;
 247 }
 248
 249
 250 // ************************ vim: set sw=4 sts=4 et: ************************ //