src/sampling/sampledSet/sampledSet/sampledSet.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 "sampledSet.H"
  28 #include "polyMesh.H"
  29 #include "primitiveMesh.H"
  30 #include "meshSearch.H"
  31 #include "writer.H"
  32
  33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37     const scalar sampledSet::tol = 1e-6;
  38
  39     defineTypeNameAndDebug(sampledSet, 0);
  40     defineRunTimeSelectionTable(sampledSet, word);
  41 }
  42
  43
  44 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  45
  46 Foam::label Foam::sampledSet::getBoundaryCell(const label faceI) const
  47 {
  48     return mesh().faceOwner()[faceI];
  49 }
  50
  51
  52 Foam::label Foam::sampledSet::getCell
  53 (
  54     const label faceI,
  55     const point& sample
  56 ) const
  57 {
  58     if (faceI == -1)
  59     {
  60         FatalErrorIn
  61         (
  62             "sampledSet::getCell(const label, const point&)"
  63         )   << "Illegal face label " << faceI
  64             << abort(FatalError);
  65     }
  66
  67     if (faceI >= mesh().nInternalFaces())
  68     {
  69         label cellI = getBoundaryCell(faceI);
  70
  71         if (!mesh().pointInCell(sample, cellI))
  72         {
  73             FatalErrorIn
  74             (
  75                 "sampledSet::getCell(const label, const point&)"
  76             )   << "Found cell " << cellI << " using face " << faceI
  77                 << ". But cell does not contain point " << sample
  78                 << abort(FatalError);
  79         }
  80         return cellI;
  81     }
  82     else
  83     {
  84         // Try owner and neighbour to see which one contains sample
  85
  86         label cellI = mesh().faceOwner()[faceI];
  87
  88         if (mesh().pointInCell(sample, cellI))
  89         {
  90             return cellI;
  91         }
  92         else
  93         {
  94             cellI = mesh().faceNeighbour()[faceI];
  95
  96             if (mesh().pointInCell(sample, cellI))
  97             {
  98                 return cellI;
  99             }
 100             else
 101             {
 102                 FatalErrorIn
 103                 (
 104                     "sampledSet::getCell(const label, const point&)"
 105                 )   << "None of the neighbours of face "
 106                     << faceI << " contains point " << sample
 107                     << abort(FatalError);
 108
 109                 return -1;
 110             }
 111         }
 112     }
 113 }
 114
 115
 116 Foam::scalar Foam::sampledSet::calcSign
 117 (
 118     const label faceI,
 119     const point& sample
 120 ) const
 121 {
 122     vector vec = sample - mesh().faceCentres()[faceI];
 123
 124     scalar magVec = mag(vec);
 125
 126     if (magVec < VSMALL)
 127     {
 128         // sample on face centre. Regard as inside
 129         return -1;
 130     }
 131
 132     vec /= magVec;
 133
 134     vector n = mesh().faceAreas()[faceI];
 135
 136     n /= mag(n) + VSMALL;
 137
 138     return n & vec;
 139 }
 140
 141
 142 // Return face (or -1) of face which is within smallDist of sample
 143 Foam::label Foam::sampledSet::findNearFace
 144 (
 145     const label cellI,
 146     const point& sample,
 147     const scalar smallDist
 148 ) const
 149 {
 150     const cell& myFaces = mesh().cells()[cellI];
 151
 152     forAll(myFaces, myFaceI)
 153     {
 154         const face& f = mesh().faces()[myFaces[myFaceI]];
 155
 156         pointHit inter = f.nearestPoint(sample, mesh().points());
 157
 158         scalar dist;
 159
 160         if (inter.hit())
 161         {
 162             dist = mag(inter.hitPoint() - sample);
 163         }
 164         else
 165         {
 166             dist = mag(inter.missPoint() - sample);
 167         }
 168
 169         if (dist < smallDist)
 170         {
 171             return myFaces[myFaceI];
 172         }
 173     }
 174     return -1;
 175 }
 176
 177
 178 // 'Pushes' point facePt (which is almost on face) in direction of cell centre
 179 // so it is clearly inside.
 180 Foam::point Foam::sampledSet::pushIn
 181 (
 182     const point& facePt,
 183     const label faceI
 184 ) const
 185 {
 186     label cellI = mesh().faceOwner()[faceI];
 187
 188     const point& cellCtr = mesh().cellCentres()[cellI];
 189
 190     point newSample =
 191         facePt + tol*(cellCtr - facePt);
 192
 193     if (!searchEngine().pointInCell(newSample, cellI))
 194     {
 195         FatalErrorIn
 196         (
 197             "sampledSet::pushIn(const point&, const label)"
 198         )   << "After pushing " << facePt << " to " << newSample
 199             << " it is still outside faceI " << faceI << endl
 200             << "Please change your starting point"
 201             << abort(FatalError);
 202     }
 203     //Info<< "pushIn : moved " << facePt << " to " << newSample
 204     //    << endl;
 205
 206     return newSample;
 207 }
 208
 209
 210 // Calculates start of tracking given samplePt and first boundary intersection
 211 // (bPoint, bFaceI). bFaceI == -1 if no boundary intersection.
 212 // Returns true if trackPt is sampling point
 213 bool Foam::sampledSet::getTrackingPoint
 214 (
 215     const vector& offset,
 216     const point& samplePt,
 217     const point& bPoint,
 218     const label bFaceI,
 219
 220     point& trackPt,
 221     label& trackCellI,
 222     label& trackFaceI
 223 ) const
 224 {
 225     const scalar smallDist = mag(tol*offset);
 226
 227     bool isGoodSample = false;
 228
 229     if (bFaceI == -1)
 230     {
 231         // No boundary intersection. Try and find cell samplePt is in
 232         trackCellI = mesh().findCell(samplePt);
 233
 234         if
 235         (
 236             (trackCellI == -1)
 237          || !mesh().pointInCell(samplePt, trackCellI)
 238         )
 239         {
 240             // Line samplePt - end_ does not intersect domain at all.
 241             // (or is along edge)
 242             //Info<< "getTrackingPoint : samplePt outside domain : "
 243             //    << "  samplePt:" << samplePt
 244             //    << endl;
 245
 246             trackCellI = -1;
 247             trackFaceI = -1;
 248
 249             isGoodSample = false;
 250         }
 251         else
 252         {
 253             // start is inside. Use it as tracking point
 254             //Info<< "getTrackingPoint : samplePt inside :"
 255             //    << "  samplePt:" << samplePt
 256             //    << "  trackCellI:" << trackCellI
 257             //    << endl;
 258
 259             trackPt = samplePt;
 260             trackFaceI = -1;
 261
 262             isGoodSample = true;
 263         }
 264     }
 265     else if (mag(samplePt - bPoint) < smallDist)
 266     {
 267         //Info<< "getTrackingPoint : samplePt:" << samplePt
 268         //    << " close to bPoint:"
 269         //    << bPoint << endl;
 270
 271         // samplePt close to bPoint. Snap to it
 272         trackPt = pushIn(bPoint, bFaceI);
 273         trackFaceI = bFaceI;
 274         trackCellI = getBoundaryCell(trackFaceI);
 275
 276         isGoodSample = true;
 277     }
 278     else
 279     {
 280         scalar sign = calcSign(bFaceI, samplePt);
 281
 282         if (sign < 0)
 283         {
 284             // samplePt inside or marginally outside.
 285             trackPt = samplePt;
 286             trackFaceI = -1;
 287             trackCellI = mesh().findCell(trackPt);
 288
 289             isGoodSample = true;
 290         }
 291         else
 292         {
 293             // samplePt outside. use bPoint
 294             trackPt = pushIn(bPoint, bFaceI);
 295             trackFaceI = bFaceI;
 296             trackCellI = getBoundaryCell(trackFaceI);
 297
 298             isGoodSample = false;
 299         }
 300     }
 301
 302     if (debug)
 303     {
 304         Info<< "sampledSet::getTrackingPoint :"
 305             << " offset:" << offset
 306             << " samplePt:" << samplePt
 307             << " bPoint:" << bPoint
 308             << " bFaceI:" << bFaceI
 309             << endl << "   Calculated first tracking point :"
 310             << " trackPt:" << trackPt
 311             << " trackCellI:" << trackCellI
 312             << " trackFaceI:" << trackFaceI
 313             << " isGoodSample:" << isGoodSample
 314             << endl;
 315     }
 316
 317     return isGoodSample;
 318 }
 319
 320
 321 void Foam::sampledSet::setSamples
 322 (
 323     const List<point>& samplingPts,
 324     const labelList& samplingCells,
 325     const labelList& samplingFaces,
 326     const labelList& samplingSegments,
 327     const scalarList& samplingCurveDist
 328 )
 329 {
 330     setSize(samplingPts.size());
 331     cells_.setSize(samplingCells.size());
 332     faces_.setSize(samplingFaces.size());
 333     segments_.setSize(samplingSegments.size());
 334     curveDist_.setSize(samplingCurveDist.size());
 335
 336     if
 337     (
 338         (cells_.size() != size())
 339      || (faces_.size() != size())
 340      || (segments_.size() != size())
 341      || (curveDist_.size() != size())
 342     )
 343     {
 344         FatalErrorIn("sampledSet::setSamples()")
 345             << "sizes not equal : "
 346             << "  points:" << size()
 347             << "  cells:" << cells_.size()
 348             << "  faces:" << faces_.size()
 349             << "  segments:" << segments_.size()
 350             << "  curveDist:" << curveDist_.size()
 351             << abort(FatalError);
 352     }
 353
 354     forAll(samplingPts, sampleI)
 355     {
 356         operator[](sampleI) = samplingPts[sampleI];
 357     }
 358     cells_ = samplingCells;
 359     faces_ = samplingFaces;
 360     segments_ = samplingSegments;
 361     curveDist_ = samplingCurveDist;
 362 }
 363
 364
 365 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 366
 367 Foam::sampledSet::sampledSet
 368 (
 369     const word& name,
 370     const polyMesh& mesh,
 371     meshSearch& searchEngine,
 372     const word& axis
 373 )
 374 :
 375     coordSet(name, axis),
 376     mesh_(mesh),
 377     searchEngine_(searchEngine),
 378     segments_(0),
 379     curveDist_(0),
 380     cells_(0),
 381     faces_(0)
 382 {}
 383
 384
 385 Foam::sampledSet::sampledSet
 386 (
 387     const word& name,
 388     const polyMesh& mesh,
 389     meshSearch& searchEngine,
 390     const dictionary& dict
 391 )
 392 :
 393     coordSet(name, dict.lookup("axis")),
 394     mesh_(mesh),
 395     searchEngine_(searchEngine),
 396     segments_(0),
 397     curveDist_(0),
 398     cells_(0),
 399     faces_(0)
 400 {}
 401
 402
 403 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 404
 405 Foam::sampledSet::~sampledSet()
 406 {}
 407
 408
 409 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 410
 411 Foam::autoPtr<Foam::sampledSet> Foam::sampledSet::New
 412 (
 413     const word& name,
 414     const polyMesh& mesh,
 415     meshSearch& searchEngine,
 416     const dictionary& dict
 417 )
 418 {
 419     word sampleType(dict.lookup("type"));
 420
 421     wordConstructorTable::iterator cstrIter =
 422         wordConstructorTablePtr_->find(sampleType);
 423
 424     if (cstrIter == wordConstructorTablePtr_->end())
 425     {
 426         FatalErrorIn
 427         (
 428             "sampledSet::New(const word&, "
 429             "const polyMesh&, meshSearch&, const dictionary&)"
 430         )   << "Unknown sample type " << sampleType
 431             << endl << endl
 432             << "Valid sample types : " << endl
 433             << wordConstructorTablePtr_->toc()
 434             << exit(FatalError);
 435     }
 436
 437     return autoPtr<sampledSet>
 438     (
 439         cstrIter()
 440         (
 441             name,
 442             mesh,
 443             searchEngine,
 444             dict
 445         )
 446     );
 447 }
 448
 449
 450 Foam::Ostream& Foam::sampledSet::write(Ostream& os) const
 451 {
 452     coordSet::write(os);
 453
 454     os  << endl << "\t(cellI)\t(faceI)" << endl;
 455
 456     forAll(*this, sampleI)
 457     {
 458         os  << '\t' << cells_[sampleI]
 459             << '\t' << faces_[sampleI]
 460             << endl;
 461     }
 462
 463     return os;
 464 }
 465
 466
 467 // ************************************************************************* //