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