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-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 "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     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     DynamicList<label>& 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 // ************************************************************************* //