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