src/meshTools/indexedOctree/treeDataPolyMeshCell.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 "treeDataPolyMeshCell.H"
  27 #include "indexedOctree.H"
  28 #include "polyMesh.H"
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 defineTypeNameAndDebug(Foam::treeDataPolyMeshCell, 0);
  33
  34
  35 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  36
  37 Foam::treeBoundBox Foam::treeDataPolyMeshCell::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::treeDataPolyMeshCell::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::treeDataPolyMeshCell::treeDataPolyMeshCell
  84 (
  85     const bool cacheBb,
  86     const polyMesh& mesh,
  87     const labelUList& cellLabels
  88 )
  89 :
  90     mesh_(mesh),
  91     cellLabels_(cellLabels),
  92     cacheBb_(cacheBb)
  93 {
  94     update();
  95 }
  96
  97
  98 Foam::treeDataPolyMeshCell::treeDataPolyMeshCell
  99 (
 100     const bool cacheBb,
 101     const polyMesh& 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::treeDataPolyMeshCell::treeDataPolyMeshCell
 114 (
 115     const bool cacheBb,
 116     const polyMesh& 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::treeDataPolyMeshCell::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::treeDataPolyMeshCell::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::treeDataPolyMeshCell::contains
 160 (
 161     const label index,
 162     const point& sample
 163 ) const
 164 {
 165     //return mesh_.pointInCell(sample, cellLabels_[index]);
 166     label tetFaceI, tetPtI;
 167     mesh_.findTetFacePt(cellLabels_[index], sample, tetFaceI, tetPtI);
 168     return tetFaceI != -1;
 169 }
 170
 171
 172 void Foam::treeDataPolyMeshCell::findNearest
 173 (
 174     const labelUList& indices,
 175     const point& sample,
 176
 177     scalar& nearestDistSqr,
 178     label& minIndex,
 179     point& nearestPoint
 180 ) const
 181 {
 182     forAll(indices, i)
 183     {
 184         label index = indices[i];
 185         label cellI = cellLabels_[index];
 186         scalar distSqr = magSqr(sample - mesh_.cellCentres()[cellI]);
 187
 188         if (distSqr < nearestDistSqr)
 189         {
 190             nearestDistSqr = distSqr;
 191             minIndex = index;
 192             nearestPoint = mesh_.cellCentres()[cellI];
 193         }
 194     }
 195 }
 196
 197
 198 bool Foam::treeDataPolyMeshCell::intersects
 199 (
 200     const label index,
 201     const point& start,
 202     const point& end,
 203     point& intersectionPoint
 204 ) const
 205 {
 206     // Do quick rejection test
 207     if (cacheBb_)
 208     {
 209         const treeBoundBox& cellBb = bbs_[index];
 210
 211         if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0)
 212         {
 213             // start and end in same block outside of cellBb.
 214             return false;
 215         }
 216     }
 217     else
 218     {
 219         const treeBoundBox cellBb = calcCellBb(cellLabels_[index]);
 220
 221         if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0)
 222         {
 223             // start and end in same block outside of cellBb.
 224             return false;
 225         }
 226     }
 227
 228
 229     // Do intersection with all faces of cell
 230     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 231
 232     // Disable picking up intersections behind us.
 233     scalar oldTol = intersection::setPlanarTol(0.0);
 234
 235     const cell& cFaces = mesh_.cells()[cellLabels_[index]];
 236
 237     const vector dir(end - start);
 238     scalar minDistSqr = magSqr(dir);
 239     bool hasMin = false;
 240
 241     forAll(cFaces, i)
 242     {
 243         const face& f = mesh_.faces()[cFaces[i]];
 244
 245         pointHit inter = f.ray
 246         (
 247             start,
 248             dir,
 249             mesh_.points(),
 250             intersection::HALF_RAY
 251         );
 252
 253         if (inter.hit() && sqr(inter.distance()) <= minDistSqr)
 254         {
 255             // Note: no extra test on whether intersection is in front of us
 256             // since using half_ray AND zero tolerance. (note that tolerance
 257             // is used to look behind us)
 258             minDistSqr = sqr(inter.distance());
 259             intersectionPoint = inter.hitPoint();
 260             hasMin = true;
 261         }
 262     }
 263
 264     // Restore picking tolerance
 265     intersection::setPlanarTol(oldTol);
 266
 267     return hasMin;
 268 }
 269
 270
 271 // ************************************************************************* //