src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.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 "ExactParticle.H"
  28
  29 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  30
  31 template<class ParticleType>
  32 template<class TrackingData>
  33 Foam::label Foam::ExactParticle<ParticleType>::track
  34 (
  35     const vector& endPosition,
  36     TrackingData& td
  37 )
  38 {
  39     this->facei_ = -1;
  40
  41     // Tracks to endPosition or stop on boundary
  42     while(!this->onBoundary() && this->stepFraction_ < 1.0 - SMALL)
  43     {
  44         this->stepFraction_ +=
  45             trackToFace(endPosition, td)
  46            *(1.0 - this->stepFraction_);
  47     }
  48
  49     return this->facei_;
  50 }
  51
  52
  53 template<class ParticleType>
  54 Foam::label Foam::ExactParticle<ParticleType>::track
  55 (
  56     const vector& endPosition
  57 )
  58 {
  59     int dummyTd;
  60     return track(endPosition, dummyTd);
  61 }
  62
  63
  64 template<class ParticleType>
  65 template<class TrackingData>
  66 Foam::scalar Foam::ExactParticle<ParticleType>::trackToFace
  67 (
  68     const vector& endPosition,
  69     TrackingData& td
  70 )
  71 {
  72     const polyMesh& mesh = this->cloud().pMesh();
  73     const labelList& cFaces = mesh.cells()[this->celli_];
  74
  75     point intersection(vector::zero);
  76     scalar trackFraction = VGREAT;
  77     label hitFacei = -1;
  78
  79     const vector vec = endPosition-this->position_;
  80
  81     forAll(cFaces, i)
  82     {
  83         label facei = cFaces[i];
  84
  85         if (facei != this->face())
  86         {
  87             pointHit inter = mesh.faces()[facei].intersection
  88             (
  89                 this->position_,
  90                 vec,
  91                 mesh.faceCentres()[facei],
  92                 mesh.points(),
  93                 intersection::HALF_RAY
  94             );
  95
  96             if (inter.hit() && inter.distance() < trackFraction)
  97             {
  98                 trackFraction = inter.distance();
  99                 hitFacei = facei;
 100                 intersection = inter.hitPoint();
 101             }
 102         }
 103     }
 104
 105     if (hitFacei == -1)
 106     {
 107         // Did not find any intersection. Fall back to original approximate
 108         // algorithm
 109         return Particle<ParticleType>::trackToFace
 110         (
 111             endPosition,
 112             td
 113         );
 114     }
 115
 116     if (trackFraction >= (1.0-SMALL))
 117     {
 118         // Nearest intersection beyond endPosition so we hit endPosition.
 119         trackFraction = 1.0;
 120         this->position_ = endPosition;
 121         this->facei_ = -1;
 122         return 1.0;
 123     }
 124     else
 125     {
 126         this->position_ = intersection;
 127         this->facei_ = hitFacei;
 128     }
 129
 130
 131     // Normal situation (trackFraction 0..1). Straight copy
 132     // of Particle::trackToFace.
 133
 134     bool internalFace = this->cloud().internalFace(this->facei_);
 135
 136     // change cell
 137     if (internalFace) // Internal face
 138     {
 139         if (this->celli_ == mesh.faceOwner()[this->facei_])
 140         {
 141             this->celli_ = mesh.faceNeighbour()[this->facei_];
 142         }
 143         else if (this->celli_ == mesh.faceNeighbour()[this->facei_])
 144         {
 145             this->celli_ = mesh.faceOwner()[this->facei_];
 146         }
 147         else
 148         {
 149             FatalErrorIn
 150             (
 151                 "ExactParticle::trackToFace"
 152                 "(const vector&, TrackingData&)"
 153             )<< "addressing failure" << nl
 154              << abort(FatalError);
 155         }
 156     }
 157     else
 158     {
 159         ParticleType& p = static_cast<ParticleType&>(*this);
 160
 161         // Soft-sphere algorithm ignores the boundary
 162         if (p.softImpact())
 163         {
 164             trackFraction = 1.0;
 165             this->position_ = endPosition;
 166         }
 167
 168         label patchi = patch(this->facei_);
 169         const polyPatch& patch = mesh.boundaryMesh()[patchi];
 170
 171         if (isA<wedgePolyPatch>(patch))
 172         {
 173             p.hitWedgePatch
 174             (
 175                 static_cast<const wedgePolyPatch&>(patch), td
 176             );
 177         }
 178         else if (isA<symmetryPolyPatch>(patch))
 179         {
 180             p.hitSymmetryPatch
 181             (
 182                 static_cast<const symmetryPolyPatch&>(patch), td
 183             );
 184         }
 185         else if (isA<cyclicPolyPatch>(patch))
 186         {
 187             p.hitCyclicPatch
 188             (
 189                 static_cast<const cyclicPolyPatch&>(patch), td
 190             );
 191         }
 192         else if (isA<processorPolyPatch>(patch))
 193         {
 194             p.hitProcessorPatch
 195             (
 196                 static_cast<const processorPolyPatch&>(patch), td
 197             );
 198         }
 199         else if (isA<wallPolyPatch>(patch))
 200         {
 201             p.hitWallPatch
 202             (
 203                 static_cast<const wallPolyPatch&>(patch), td
 204             );
 205         }
 206         else if (isA<polyPatch>(patch))
 207         {
 208             p.hitPatch
 209             (
 210                 static_cast<const polyPatch&>(patch), td
 211             );
 212         }
 213         else
 214         {
 215             FatalErrorIn
 216             (
 217                 "ExactParticle::trackToFace"
 218                 "(const vector& endPosition, scalar& trackFraction)"
 219             )<< "patch type " << patch.type() << " not suported" << nl
 220              << abort(FatalError);
 221         }
 222     }
 223
 224     // If the trackFraction = 0 something went wrong.
 225     // Either the particle is flipping back and forth across a face perhaps
 226     // due to velocity interpolation errors or it is in a "hole" in the mesh
 227     // caused by face warpage.
 228     // In both cases resolve the positional ambiguity by moving the particle
 229     // slightly towards the cell-centre.
 230     if (trackFraction < SMALL)
 231     {
 232         this->position_ +=
 233             1.0e-6*(mesh.cellCentres()[this->celli_] - this->position_);
 234     }
 235
 236     return trackFraction;
 237 }
 238
 239
 240 template<class ParticleType>
 241 Foam::scalar Foam::ExactParticle<ParticleType>::trackToFace
 242 (
 243     const vector& endPosition
 244 )
 245 {
 246     int dummyTd;
 247     return trackToFace(endPosition, dummyTd);
 248 }
 249
 250
 251 template<class ParticleType>
 252 Foam::Ostream& Foam::operator<<
 253 (
 254     Ostream& os,
 255     const ExactParticle<ParticleType>& p
 256 )
 257 {
 258     return operator<<(os, static_cast<const Particle<ParticleType>&>(p));
 259 }
 260
 261
 262 // ************************************************************************* //