src/meshTools/searchableSurface/searchableSphere.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 "searchableSphere.H"
  27 #include "addToRunTimeSelectionTable.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33
  34 defineTypeNameAndDebug(searchableSphere, 0);
  35 addToRunTimeSelectionTable(searchableSurface, searchableSphere, dict);
  36
  37 }
  38
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 Foam::pointIndexHit Foam::searchableSphere::findNearest
  43 (
  44     const point& sample,
  45     const scalar nearestDistSqr
  46 ) const
  47 {
  48     pointIndexHit info(false, sample, -1);
  49
  50     const vector n(sample-centre_);
  51     scalar magN = mag(n);
  52
  53     if (nearestDistSqr > sqr(magN-radius_))
  54     {
  55         if (magN < ROOTVSMALL)
  56         {
  57             info.rawPoint() = centre_ + vector(1,0,0)/magN*radius_;
  58         }
  59         else
  60         {
  61             info.rawPoint() = centre_ + n/magN*radius_;
  62         }
  63         info.setHit();
  64         info.setIndex(0);
  65     }
  66
  67     return info;
  68 }
  69
  70
  71 // From Graphics Gems - intersection of sphere with ray
  72 void Foam::searchableSphere::findLineAll
  73 (
  74     const point& start,
  75     const point& end,
  76     pointIndexHit& near,
  77     pointIndexHit& far
  78 ) const
  79 {
  80     near.setMiss();
  81     far.setMiss();
  82
  83     vector dir(end-start);
  84     scalar magSqrDir = magSqr(dir);
  85
  86     if (magSqrDir > ROOTVSMALL)
  87     {
  88         const vector toCentre(centre_-start);
  89         scalar magSqrToCentre = magSqr(toCentre);
  90
  91         dir /= Foam::sqrt(magSqrDir);
  92
  93         scalar v = (toCentre & dir);
  94
  95         scalar disc = sqr(radius_) - (magSqrToCentre - sqr(v));
  96
  97         if (disc >= 0)
  98         {
  99             scalar d = Foam::sqrt(disc);
 100
 101             scalar nearParam = v-d;
 102
 103             if (nearParam >= 0 && sqr(nearParam) <= magSqrDir)
 104             {
 105                 near.setHit();
 106                 near.setPoint(start + nearParam*dir);
 107                 near.setIndex(0);
 108             }
 109
 110             scalar farParam = v+d;
 111
 112             if (farParam >= 0 && sqr(farParam) <= magSqrDir)
 113             {
 114                 far.setHit();
 115                 far.setPoint(start + farParam*dir);
 116                 far.setIndex(0);
 117             }
 118         }
 119     }
 120 }
 121
 122
 123 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 124
 125 Foam::searchableSphere::searchableSphere
 126 (
 127     const IOobject& io,
 128     const point& centre,
 129     const scalar radius
 130 )
 131 :
 132     searchableSurface(io),
 133     centre_(centre),
 134     radius_(radius)
 135 {}
 136
 137
 138 Foam::searchableSphere::searchableSphere
 139 (
 140     const IOobject& io,
 141     const dictionary& dict
 142 )
 143 :
 144     searchableSurface(io),
 145     centre_(dict.lookup("centre")),
 146     radius_(readScalar(dict.lookup("radius")))
 147 {}
 148
 149
 150 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 151
 152 Foam::searchableSphere::~searchableSphere()
 153 {}
 154
 155
 156 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 157
 158 const Foam::wordList& Foam::searchableSphere::regions() const
 159 {
 160     if (regions_.empty())
 161     {
 162         regions_.setSize(1);
 163         regions_[0] = "region0";
 164     }
 165     return regions_;
 166 }
 167
 168
 169
 170 void Foam::searchableSphere::findNearest
 171 (
 172     const pointField& samples,
 173     const scalarField& nearestDistSqr,
 174     List<pointIndexHit>& info
 175 ) const
 176 {
 177     info.setSize(samples.size());
 178
 179     forAll(samples, i)
 180     {
 181         info[i] = findNearest(samples[i], nearestDistSqr[i]);
 182     }
 183 }
 184
 185
 186 void Foam::searchableSphere::findLine
 187 (
 188     const pointField& start,
 189     const pointField& end,
 190     List<pointIndexHit>& info
 191 ) const
 192 {
 193     info.setSize(start.size());
 194
 195     pointIndexHit b;
 196
 197     forAll(start, i)
 198     {
 199         // Pick nearest intersection. If none intersected take second one.
 200         findLineAll(start[i], end[i], info[i], b);
 201         if (!info[i].hit() && b.hit())
 202         {
 203             info[i] = b;
 204         }
 205     }
 206 }
 207
 208
 209 void Foam::searchableSphere::findLineAny
 210 (
 211     const pointField& start,
 212     const pointField& end,
 213     List<pointIndexHit>& info
 214 ) const
 215 {
 216     info.setSize(start.size());
 217
 218     pointIndexHit b;
 219
 220     forAll(start, i)
 221     {
 222         // Discard far intersection
 223         findLineAll(start[i], end[i], info[i], b);
 224         if (!info[i].hit() && b.hit())
 225         {
 226             info[i] = b;
 227         }
 228     }
 229 }
 230
 231
 232 void Foam::searchableSphere::findLineAll
 233 (
 234     const pointField& start,
 235     const pointField& end,
 236     List<List<pointIndexHit> >& info
 237 ) const
 238 {
 239     info.setSize(start.size());
 240
 241     forAll(start, i)
 242     {
 243         pointIndexHit near, far;
 244         findLineAll(start[i], end[i], near, far);
 245
 246         if (near.hit())
 247         {
 248             if (far.hit())
 249             {
 250                 info[i].setSize(2);
 251                 info[i][0] = near;
 252                 info[i][1] = far;
 253             }
 254             else
 255             {
 256                 info[i].setSize(1);
 257                 info[i][0] = near;
 258             }
 259         }
 260         else
 261         {
 262             if (far.hit())
 263             {
 264                 info[i].setSize(1);
 265                 info[i][0] = far;
 266             }
 267             else
 268             {
 269                 info[i].clear();
 270             }
 271         }
 272     }
 273 }
 274
 275
 276 void Foam::searchableSphere::getRegion
 277 (
 278     const List<pointIndexHit>& info,
 279     labelList& region
 280 ) const
 281 {
 282     region.setSize(info.size());
 283     region = 0;
 284 }
 285
 286
 287 void Foam::searchableSphere::getNormal
 288 (
 289     const List<pointIndexHit>& info,
 290     vectorField& normal
 291 ) const
 292 {
 293     normal.setSize(info.size());
 294     normal = vector::zero;
 295
 296     forAll(info, i)
 297     {
 298         if (info[i].hit())
 299         {
 300             normal[i] = info[i].hitPoint() - centre_;
 301
 302             normal[i] /= mag(normal[i])+VSMALL;
 303         }
 304         else
 305         {
 306             // Set to what?
 307         }
 308     }
 309 }
 310
 311
 312 void Foam::searchableSphere::getVolumeType
 313 (
 314     const pointField& points,
 315     List<volumeType>& volType
 316 ) const
 317 {
 318     volType.setSize(points.size());
 319     volType = INSIDE;
 320
 321     forAll(points, pointI)
 322     {
 323         const point& pt = points[pointI];
 324
 325         if (magSqr(pt - centre_) <= sqr(radius_))
 326         {
 327             volType[pointI] = INSIDE;
 328         }
 329         else
 330         {
 331             volType[pointI] = OUTSIDE;
 332         }
 333     }
 334 }
 335
 336
 337 // ************************************************************************* //