src/lagrangian/basic/Particle/Particle.C

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