src/lagrangian/basic/Particle/Particle.C

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