src/lagrangian/basic/Particle/Particle.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 "Particle.H"
  28 #include "Cloud.H"
  29 #include "wedgePolyPatch.H"
  30 #include "symmetryPolyPatch.H"
  31 #include "cyclicPolyPatch.H"
  32 #include "processorPolyPatch.H"
  33 #include "wallPolyPatch.H"
  34 #include "transform.H"
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38 template<class ParticleType>
  39 Foam::labelList Foam::Particle<ParticleType>::findFaces
  40 (
  41     const vector& position
  42 ) const
  43 {
  44     const polyMesh& mesh = cloud_.polyMesh_;
  45     const labelList& faces = mesh.cells()[celli_];
  46     const vector& C = mesh.cellCentres()[celli_];
  47
  48     labelList faceList(0);
  49     forAll(faces, i)
  50     {
  51         label facei = faces[i];
  52         scalar lam = lambda(C, position, facei);
  53
  54         if ((lam > 0) && (lam < 1.0))
  55         {
  56             label n = faceList.size();
  57             faceList.setSize(n+1);
  58             faceList[n] = facei;
  59         }
  60     }
  61
  62     return faceList;
  63 }
  64
  65
  66 template<class ParticleType>
  67 Foam::labelList Foam::Particle<ParticleType>::findFaces
  68 (
  69     const vector& position,
  70     const label celli,
  71     const scalar stepFraction
  72 ) const
  73 {
  74     const polyMesh& mesh = cloud_.polyMesh_;
  75     const labelList& faces = mesh.cells()[celli];
  76     const vector& C = mesh.cellCentres()[celli];
  77
  78     labelList faceList(0);
  79     forAll(faces, i)
  80     {
  81         label facei = faces[i];
  82         scalar lam = lambda(C, position, facei, stepFraction);
  83
  84         if ((lam > 0) && (lam < 1.0))
  85         {
  86             label n = faceList.size();
  87             faceList.setSize(n+1);
  88             faceList[n] = facei;
  89         }
  90     }
  91
  92     return faceList;
  93 }
  94
  95
  96 template<class ParticleType>
  97 template<class TrackData>
  98 void Foam::Particle<ParticleType>::prepareForParallelTransfer
  99 (
 100     const label patchi,
 101     TrackData& td
 102 )
 103 {
 104     // Convert the face index to be local to the processor patch
 105     facei_ = patchFace(patchi, facei_);
 106 }
 107
 108
 109 template<class ParticleType>
 110 template<class TrackData>
 111 void Foam::Particle<ParticleType>::correctAfterParallelTransfer
 112 (
 113     const label patchi,
 114     TrackData& td
 115 )
 116 {
 117     const processorPolyPatch& ppp =
 118         refCast<const processorPolyPatch>
 119         (cloud_.pMesh().boundaryMesh()[patchi]);
 120
 121     celli_ = ppp.faceCells()[facei_];
 122
 123     if (!ppp.parallel())
 124     {
 125         if (ppp.forwardT().size() == 1)
 126         {
 127             const tensor& T = ppp.forwardT()[0];
 128             transformPosition(T);
 129             static_cast<ParticleType&>(*this).transformProperties(T);
 130         }
 131         else
 132         {
 133             const tensor& T = ppp.forwardT()[facei_];
 134             transformPosition(T);
 135             static_cast<ParticleType&>(*this).transformProperties(T);
 136         }
 137     }
 138     else if (ppp.separated())
 139     {
 140         if (ppp.separation().size() == 1)
 141         {
 142             position_ -= ppp.separation()[0];
 143             static_cast<ParticleType&>(*this).transformProperties
 144             (
 145                 -ppp.separation()[0]
 146             );
 147         }
 148         else
 149         {
 150             position_ -= ppp.separation()[facei_];
 151             static_cast<ParticleType&>(*this).transformProperties
 152             (
 153                 -ppp.separation()[facei_]
 154             );
 155         }
 156     }
 157
 158     // Reset the face index for the next tracking operation
 159     if (stepFraction_ > (1.0 - SMALL))
 160     {
 161         stepFraction_ = 1.0;
 162         facei_ = -1;
 163     }
 164     else
 165     {
 166         facei_ += ppp.start();
 167     }
 168 }
 169
 170
 171 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 172
 173 template<class ParticleType>
 174 Foam::Particle<ParticleType>::Particle
 175 (
 176     const Cloud<ParticleType>& cloud,
 177     const vector& position,
 178     const label celli
 179 )
 180 :
 181     cloud_(cloud),
 182     position_(position),
 183     celli_(celli),
 184     facei_(-1),
 185     stepFraction_(0.0)
 186 {}
 187
 188
 189 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 190
 191 template<class ParticleType>
 192 template<class TrackData>
 193 Foam::label Foam::Particle<ParticleType>::track
 194 (
 195     const vector& endPosition,
 196     TrackData& td
 197 )
 198 {
 199     facei_ = -1;
 200
 201     // Tracks to endPosition or stop on boundary
 202     while(!onBoundary() && stepFraction_ < 1.0 - SMALL)
 203     {
 204         stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);
 205     }
 206
 207     return facei_;
 208 }
 209
 210
 211
 212 template<class ParticleType>
 213 Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
 214 {
 215     int dummyTd;
 216     return track(endPosition, dummyTd);
 217 }
 218
 219 template<class ParticleType>
 220 template<class TrackData>
 221 Foam::scalar Foam::Particle<ParticleType>::trackToFace
 222 (
 223     const vector& endPosition,
 224     TrackData& td
 225 )
 226 {
 227     const polyMesh& mesh = cloud_.polyMesh_;
 228
 229     labelList faces = findFaces(endPosition);
 230
 231     facei_ = -1;
 232     scalar trackFraction = 0.0;
 233
 234     if (faces.size() == 0) // inside cell
 235     {
 236         trackFraction = 1.0;
 237         position_ = endPosition;
 238     }
 239     else // hit face
 240     {
 241         scalar lambdaMin = GREAT;
 242
 243         if (faces.size() == 1)
 244         {
 245             lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
 246             facei_ = faces[0];
 247         }
 248         else
 249         {
 250             // If the particle has to cross more than one cell to reach the
 251             // endPosition, we check which way to go.
 252             // If one of the faces is a boundary face and the particle is
 253             // outside, we choose the boundary face.
 254             // The particle is outside if one of the lambda's is > 1 or < 0
 255             forAll(faces, i)
 256             {
 257                 scalar lam =
 258                     lambda(position_, endPosition, faces[i], stepFraction_);
 259
 260                 if (lam < lambdaMin)
 261                 {
 262                     lambdaMin = lam;
 263                     facei_ = faces[i];
 264                 }
 265             }
 266         }
 267
 268         bool internalFace = cloud_.internalFace(facei_);
 269
 270         // For warped faces the particle can be 'outside' the cell.
 271         // This will yield a lambda larger than 1, or smaller than 0
 272         // For values < 0, the particle travels away from the cell
 273         // and we don't move the particle, only change cell.
 274         // For values larger than 1, we move the particle to endPosition only.
 275         if (lambdaMin > 0.0)
 276         {
 277             if (lambdaMin <= 1.0)
 278             {
 279                 trackFraction = lambdaMin;
 280                 position_ += trackFraction*(endPosition - position_);
 281             }
 282             else
 283             {
 284                 trackFraction = 1.0;
 285                 position_ = endPosition;
 286             }
 287         }
 288         else if (static_cast<ParticleType&>(*this).softImpact())
 289         {
 290             // Soft-sphere particles can travel outside the domain
 291             // but we don't use lambda since this the particle
 292             // is going away from face
 293             trackFraction = 1.0;
 294             position_ = endPosition;
 295         }
 296
 297         // change cell
 298         if (internalFace) // Internal face
 299         {
 300             if (celli_ == cloud_.owner_[facei_])
 301             {
 302                 celli_ = cloud_.neighbour_[facei_];
 303             }
 304             else if (celli_ == cloud_.neighbour_[facei_])
 305             {
 306                 celli_ = cloud_.owner_[facei_];
 307             }
 308             else
 309             {
 310                 FatalErrorIn
 311                 (
 312                     "Particle::trackToFace(const vector&, TrackData&)"
 313                 )<< "addressing failure" << nl
 314                  << abort(FatalError);
 315             }
 316         }
 317         else
 318         {
 319             ParticleType& p = static_cast<ParticleType&>(*this);
 320
 321             // Soft-sphere algorithm ignores the boundary
 322             if (p.softImpact())
 323             {
 324                 trackFraction = 1.0;
 325                 position_ = endPosition;
 326             }
 327
 328             label patchi = patch(facei_);
 329             const polyPatch& patch = mesh.boundaryMesh()[patchi];
 330
 331             if (isA<wedgePolyPatch>(patch))
 332             {
 333                 p.hitWedgePatch
 334                 (
 335                     static_cast<const wedgePolyPatch&>(patch), td
 336                 );
 337             }
 338             else if (isA<symmetryPolyPatch>(patch))
 339             {
 340                 p.hitSymmetryPatch
 341                 (
 342                     static_cast<const symmetryPolyPatch&>(patch), td
 343                 );
 344             }
 345             else if (isA<cyclicPolyPatch>(patch))
 346             {
 347                 p.hitCyclicPatch
 348                 (
 349                     static_cast<const cyclicPolyPatch&>(patch), td
 350                 );
 351             }
 352             else if (isA<processorPolyPatch>(patch))
 353             {
 354                 p.hitProcessorPatch
 355                 (
 356                     static_cast<const processorPolyPatch&>(patch), td
 357                 );
 358             }
 359             else if (isA<wallPolyPatch>(patch))
 360             {
 361                 p.hitWallPatch
 362                 (
 363                     static_cast<const wallPolyPatch&>(patch), td
 364                 );
 365             }
 366             else if (isA<polyPatch>(patch))
 367             {
 368                 p.hitPatch
 369                 (
 370                     static_cast<const polyPatch&>(patch), td
 371                 );
 372             }
 373             else
 374             {
 375                 FatalErrorIn
 376                 (
 377                     "Particle::trackToFace"
 378                     "(const vector& endPosition, scalar& trackFraction)"
 379                 )<< "patch type " << patch.type() << " not suported" << nl
 380                  << abort(FatalError);
 381             }
 382         }
 383     }
 384
 385     // If the trackFraction = 0 something went wrong.
 386     // Either the particle is flipping back and forth across a face perhaps
 387     // due to velocity interpolation errors or it is in a "hole" in the mesh
 388     // caused by face warpage.
 389     // In both cases resolve the positional ambiguity by moving the particle
 390     // slightly towards the cell-centre.
 391     if (trackFraction < SMALL)
 392     {
 393         position_ += 1.0e-6*(mesh.cellCentres()[celli_] - position_);
 394     }
 395
 396     return trackFraction;
 397 }
 398
 399 template<class ParticleType>
 400 Foam::scalar Foam::Particle<ParticleType>::trackToFace
 401 (
 402     const vector& endPosition
 403 )
 404 {
 405     int dummyTd;
 406     return trackToFace(endPosition, dummyTd);
 407 }
 408
 409 template<class ParticleType>
 410 void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
 411 {
 412     position_ = transform(T, position_);
 413 }
 414
 415
 416 template<class ParticleType>
 417 void Foam::Particle<ParticleType>::transformProperties(const tensor&)
 418 {}
 419
 420
 421 template<class ParticleType>
 422 void Foam::Particle<ParticleType>::transformProperties(const vector&)
 423 {}
 424
 425
 426 template<class ParticleType>
 427 template<class TrackData>
 428 void Foam::Particle<ParticleType>::hitWedgePatch
 429 (
 430     const wedgePolyPatch& wpp,
 431     TrackData&
 432 )
 433 {
 434     vector nf = wpp.faceAreas()[wpp.whichFace(facei_)];
 435     nf /= mag(nf);
 436
 437     static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
 438 }
 439
 440
 441 template<class ParticleType>
 442 template<class TrackData>
 443 void Foam::Particle<ParticleType>::hitSymmetryPatch
 444 (
 445     const symmetryPolyPatch& spp,
 446     TrackData&
 447 )
 448 {
 449     vector nf = spp.faceAreas()[spp.whichFace(facei_)];
 450     nf /= mag(nf);
 451
 452     static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
 453 }
 454
 455
 456 template<class ParticleType>
 457 template<class TrackData>
 458 void Foam::Particle<ParticleType>::hitCyclicPatch
 459 (
 460     const cyclicPolyPatch& cpp,
 461     TrackData&
 462 )
 463 {
 464     label patchFacei_ = cpp.whichFace(facei_);
 465
 466     facei_ = cpp.transformGlobalFace(facei_);
 467
 468     celli_ = cloud_.polyMesh_.faceOwner()[facei_];
 469
 470     if (!cpp.parallel())
 471     {
 472         const tensor& T = cpp.transformT(patchFacei_);
 473
 474         transformPosition(T);
 475         static_cast<ParticleType&>(*this).transformProperties(T);
 476     }
 477     else if (cpp.separated())
 478     {
 479         position_ += cpp.separation(patchFacei_);
 480         static_cast<ParticleType&>(*this).transformProperties
 481         (
 482             cpp.separation(patchFacei_)
 483         );
 484     }
 485 }
 486
 487
 488 template<class ParticleType>
 489 template<class TrackData>
 490 void Foam::Particle<ParticleType>::hitProcessorPatch
 491 (
 492     const processorPolyPatch& spp,
 493     TrackData& td
 494 )
 495 {}
 496
 497
 498 template<class ParticleType>
 499 template<class TrackData>
 500 void Foam::Particle<ParticleType>::hitWallPatch
 501 (
 502     const wallPolyPatch& spp,
 503     TrackData&
 504 )
 505 {}
 506
 507
 508 template<class ParticleType>
 509 template<class TrackData>
 510 void Foam::Particle<ParticleType>::hitPatch
 511 (
 512     const polyPatch& spp,
 513     TrackData&
 514 )
 515 {}
 516
 517
 518 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 519
 520 #include "ParticleIO.C"
 521
 522 // ************************************************************************* //