src/lagrangian/basic/particle/particleI.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "polyMesh.H"
  27 #include "Time.H"
  28
  29 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  30
  31 inline void Foam::particle::findTris
  32 (
  33     const vector& position,
  34     DynamicList<label>& faceList,
  35     const tetPointRef& tet,
  36     const FixedList<vector, 4>& tetAreas,
  37     const FixedList<label, 4>& tetPlaneBasePtIs,
  38     const scalar tol
  39 ) const
  40 {
  41     faceList.clear();
  42
  43     const point Ct = tet.centre();
  44
  45     for (label i = 0; i < 4; i++)
  46     {
  47         scalar lambda = tetLambda
  48         (
  49             Ct,
  50             position,
  51             i,
  52             tetAreas[i],
  53             tetPlaneBasePtIs[i],
  54             cellI_,
  55             tetFaceI_,
  56             tetPtI_,
  57             tol
  58         );
  59
  60         if ((lambda > 0.0) && (lambda < 1.0))
  61         {
  62             faceList.append(i);
  63         }
  64     }
  65 }
  66
  67
  68 inline Foam::scalar Foam::particle::tetLambda
  69 (
  70     const vector& from,
  71     const vector& to,
  72     const label triI,
  73     const vector& n,
  74     const label tetPlaneBasePtI,
  75     const label cellI,
  76     const label tetFaceI,
  77     const label tetPtI,
  78     const scalar tol
  79 ) const
  80 {
  81     const pointField& pPts = mesh_.points();
  82
  83     if (mesh_.moving())
  84     {
  85         return movingTetLambda
  86         (
  87             from,
  88             to,
  89             triI,
  90             n,
  91             tetPlaneBasePtI,
  92             cellI,
  93             tetFaceI,
  94             tetPtI,
  95             tol
  96         );
  97     }
  98
  99     const point& base = pPts[tetPlaneBasePtI];
 100
 101     scalar lambdaNumerator = (base - from) & n;
 102     scalar lambdaDenominator = (to - from) & n;
 103
 104     // n carries the area of the tet faces, so the dot product with a
 105     // delta-length has the units of volume.  Comparing the component of each
 106     // delta-length in the direction of n times the face area to a fraction of
 107     // the cell volume.
 108
 109     if (mag(lambdaDenominator) < tol)
 110     {
 111         if (mag(lambdaNumerator) < tol)
 112         {
 113             // Track starts on the face, and is potentially
 114             // parallel to it.  +-tol/+-tol is not a good
 115             // comparison, return 0.0, in anticipation of tet
 116             // centre correction.
 117
 118             return 0.0;
 119         }
 120         else
 121         {
 122             if (mag((to - from)) < tol/mag(n))
 123             {
 124                 // 'Zero' length track (compared to the tolerance, which is
 125                 // based on the cell volume, divided by the tet face area), not
 126                 // along the face, face cannot be crossed.
 127
 128                 return GREAT;
 129             }
 130             else
 131             {
 132                 // Trajectory is non-zero and parallel to face
 133
 134                 lambdaDenominator = sign(lambdaDenominator)*SMALL;
 135             }
 136         }
 137     }
 138
 139     return lambdaNumerator/lambdaDenominator;
 140 }
 141
 142
 143 inline Foam::scalar Foam::particle::movingTetLambda
 144 (
 145     const vector& from,
 146     const vector& to,
 147     const label triI,
 148     const vector& n,
 149     const label tetPlaneBasePtI,
 150     const label cellI,
 151     const label tetFaceI,
 152     const label tetPtI,
 153     const scalar tol
 154 ) const
 155 {
 156     const pointField& pPts = mesh_.points();
 157     const pointField& oldPPts = mesh_.oldPoints();
 158
 159     // Base point of plane at end of motion
 160     const point& b = pPts[tetPlaneBasePtI];
 161
 162     // n: Normal of plane at end of motion
 163
 164     // Base point of plane at start of timestep
 165     const point& b00 = oldPPts[tetPlaneBasePtI];
 166
 167     // Base point of plane at start of tracking portion (cast forward by
 168     // stepFraction)
 169     point b0 = b00 + stepFraction_*(b - b00);
 170
 171     // Normal of plane at start of tracking portion
 172     vector n0 = vector::zero;
 173
 174     {
 175         tetIndices tetIs(cellI, tetFaceI, tetPtI, mesh_);
 176
 177         // tet at timestep start
 178         tetPointRef tet00 = tetIs.oldTet(mesh_);
 179
 180         // tet at timestep end
 181         tetPointRef tet = tetIs.tet(mesh_);
 182
 183         point tet0PtA = tet00.a() + stepFraction_*(tet.a() - tet00.a());
 184         point tet0PtB = tet00.b() + stepFraction_*(tet.b() - tet00.b());
 185         point tet0PtC = tet00.c() + stepFraction_*(tet.c() - tet00.c());
 186         point tet0PtD = tet00.d() + stepFraction_*(tet.d() - tet00.d());
 187
 188         // Tracking portion start tet (cast forward by stepFraction)
 189         tetPointRef tet0(tet0PtA, tet0PtB, tet0PtC, tet0PtD);
 190
 191         switch (triI)
 192         {
 193             case 0:
 194             {
 195                 n0 = tet0.Sa();
 196                 break;
 197             }
 198             case 1:
 199             {
 200                 n0 = tet0.Sb();
 201                 break;
 202             }
 203             case 2:
 204             {
 205                 n0 = tet0.Sc();
 206                 break;
 207             }
 208             case 3:
 209             {
 210                 n0 = tet0.Sd();
 211                 break;
 212             }
 213             default:
 214             {
 215                 break;
 216             }
 217         }
 218     }
 219
 220     if (mag(n0) < SMALL)
 221     {
 222         // If the old normal is zero (for example in layer addition)
 223         // then use the current normal;
 224
 225         n0 = n;
 226     }
 227
 228     scalar lambdaNumerator = 0;
 229     scalar lambdaDenominator = 0;
 230
 231     vector dP = to - from;
 232     vector dN = n - n0;
 233     vector dB = b - b0;
 234     vector dS = from - b0;
 235
 236     if (mag(dN) > SMALL)
 237     {
 238         scalar a = (dP - dB) & dN;
 239         scalar b = ((dP - dB) & n0) + (dS & dN);
 240         scalar c = dS & n0;
 241
 242         if (mag(a) > SMALL)
 243         {
 244
 245             // Solve quadratic for lambda
 246             scalar discriminant = sqr(b) - 4.0*a*c;
 247
 248             if (discriminant < 0)
 249             {
 250                 // Imaginary roots only - face not crossed
 251                 return GREAT;
 252             }
 253             else
 254             {
 255                 // Numerical Recipes in C,
 256                 // Second Edition (1992),
 257                 // Section 5.6.
 258                 // q = -0.5*(b + sgn(b)*sqrt(b^2 - 4ac))
 259                 // x1 = q/a
 260                 // x2 = c/q
 261
 262                 scalar q = -0.5*(b + sign(b)*Foam::sqrt(discriminant));
 263
 264                 if (mag(q) < VSMALL)
 265                 {
 266                     // If q is zero, then l1 = q/a is the required
 267                     // value of lambda, and is zero.
 268
 269                     return 0.0;
 270                 }
 271
 272                 scalar l1 = q/a;
 273                 scalar l2 = c/q;
 274
 275                 // There will be two roots, a big one and a little
 276                 // one, choose the little one.
 277
 278                 if (mag(l1) < mag(l2))
 279                 {
 280                     return l1;
 281                 }
 282                 else
 283                 {
 284                     return l2;
 285                 }
 286             }
 287         }
 288         {
 289             // When a is zero, solve the first order polynomial
 290
 291             lambdaNumerator = -c;
 292             lambdaDenominator = b;
 293         }
 294     }
 295     else
 296     {
 297         // when n = n0 is zero, there is no plane rotation, solve the
 298         // first order polynomial
 299
 300         lambdaNumerator = -(dS & n0);
 301         lambdaDenominator = ((dP - dB) & n0);
 302
 303     }
 304
 305     if (mag(lambdaDenominator) < tol)
 306     {
 307         if (mag(lambdaNumerator) < tol)
 308         {
 309             // Track starts on the face, and is potentially
 310             // parallel to it.  +-tol)/+-tol is not a good
 311             // comparison, return 0.0, in anticipation of tet
 312             // centre correction.
 313
 314             return 0.0;
 315         }
 316         else
 317         {
 318             if (mag((to - from)) < tol/mag(n))
 319             {
 320                 // Zero length track, not along the face, face
 321                 // cannot be crossed.
 322
 323                 return GREAT;
 324             }
 325             else
 326             {
 327                 // Trajectory is non-zero and parallel to face
 328
 329                 lambdaDenominator = sign(lambdaDenominator)*SMALL;
 330             }
 331         }
 332     }
 333
 334     return lambdaNumerator/lambdaDenominator;
 335 }
 336
 337
 338
 339 inline void Foam::particle::tetNeighbour(label triI)
 340 {
 341     const labelList& pOwner = mesh_.faceOwner();
 342     const faceList& pFaces = mesh_.faces();
 343
 344     bool own = (pOwner[tetFaceI_] == cellI_);
 345
 346     const Foam::face& f = pFaces[tetFaceI_];
 347
 348     label tetBasePtI = mesh_.tetBasePtIs()[tetFaceI_];
 349
 350     if (tetBasePtI == -1)
 351     {
 352         FatalErrorIn
 353         (
 354             "inline void Foam::particle::tetNeighbour(label triI)"
 355         )
 356             << "No base point for face " << tetFaceI_ << ", " << f
 357             << ", produces a valid tet decomposition."
 358             << abort(FatalError);
 359     }
 360
 361     label facePtI = (tetPtI_ + tetBasePtI) % f.size();
 362     label otherFacePtI = f.fcIndex(facePtI);
 363
 364     switch (triI)
 365     {
 366         case 0:
 367         {
 368             // Crossing this triangle changes tet to that in the
 369             // neighbour cell over tetFaceI
 370
 371             // Modification of cellI_ will happen by other indexing,
 372             // tetFaceI_ and tetPtI don't change.
 373
 374             break;
 375         }
 376         case 1:
 377         {
 378             crossEdgeConnectedFace
 379             (
 380                 cellI_,
 381                 tetFaceI_,
 382                 tetPtI_,
 383                 Foam::edge(f[facePtI], f[otherFacePtI])
 384             );
 385
 386             break;
 387         }
 388         case 2:
 389         {
 390             if (own)
 391             {
 392                 if (tetPtI_ < f.size() - 2)
 393                 {
 394                     tetPtI_ = f.fcIndex(tetPtI_);
 395                 }
 396                 else
 397                 {
 398                     crossEdgeConnectedFace
 399                     (
 400                         cellI_,
 401                         tetFaceI_,
 402                         tetPtI_,
 403                         Foam::edge(f[tetBasePtI], f[otherFacePtI])
 404                     );
 405                 }
 406             }
 407             else
 408             {
 409                 if (tetPtI_ > 1)
 410                 {
 411                     tetPtI_ = f.rcIndex(tetPtI_);
 412                 }
 413                 else
 414                 {
 415                     crossEdgeConnectedFace
 416                     (
 417                         cellI_,
 418                         tetFaceI_,
 419                         tetPtI_,
 420                         Foam::edge(f[tetBasePtI], f[facePtI])
 421                     );
 422                 }
 423             }
 424
 425             break;
 426         }
 427         case 3:
 428         {
 429             if (own)
 430             {
 431                 if (tetPtI_ > 1)
 432                 {
 433                     tetPtI_ = f.rcIndex(tetPtI_);
 434                 }
 435                 else
 436                 {
 437                     crossEdgeConnectedFace
 438                     (
 439                         cellI_,
 440                         tetFaceI_,
 441                         tetPtI_,
 442                         Foam::edge(f[tetBasePtI], f[facePtI])
 443                     );
 444                 }
 445             }
 446             else
 447             {
 448                 if (tetPtI_ < f.size() - 2)
 449                 {
 450                     tetPtI_ = f.fcIndex(tetPtI_);
 451                 }
 452                 else
 453                 {
 454                     crossEdgeConnectedFace
 455                     (
 456                         cellI_,
 457                         tetFaceI_,
 458                         tetPtI_,
 459                         Foam::edge(f[tetBasePtI], f[otherFacePtI])
 460                     );
 461                 }
 462             }
 463
 464             break;
 465         }
 466         default:
 467         {
 468             FatalErrorIn
 469             (
 470                 "inline void "
 471                 "Foam::particle::tetNeighbour(label triI)"
 472             )
 473                 << "Tet tri face index error, can only be 0..3, supplied "
 474                 << triI << abort(FatalError);
 475
 476             break;
 477         }
 478     }
 479 }
 480
 481
 482 inline void Foam::particle::crossEdgeConnectedFace
 483 (
 484     const label& cellI,
 485     label& tetFaceI,
 486     label& tetPtI,
 487     const edge& e
 488 )
 489 {
 490     const faceList& pFaces = mesh_.faces();
 491     const cellList& pCells = mesh_.cells();
 492
 493     const Foam::face& f = pFaces[tetFaceI];
 494
 495     const Foam::cell& thisCell = pCells[cellI];
 496
 497     forAll(thisCell, cFI)
 498     {
 499         // Loop over all other faces of this cell and
 500         // find the one that shares this edge
 501
 502         label fI = thisCell[cFI];
 503
 504         if (tetFaceI == fI)
 505         {
 506             continue;
 507         }
 508
 509         const Foam::face& otherFace = pFaces[fI];
 510
 511         label edDir = otherFace.edgeDirection(e);
 512
 513         if (edDir == 0)
 514         {
 515             continue;
 516         }
 517         else if (f == pFaces[fI])
 518         {
 519             // This is a necessary condition if using duplicate baffles
 520             // (so coincident faces). We need to make sure we don't cross into
 521             // the face with the same vertices since we might enter a tracking
 522             // loop where it never exits. This test should be cheap
 523             // for most meshes so can be left in for 'normal' meshes.
 524
 525             continue;
 526         }
 527         else
 528         {
 529             //Found edge on other face
 530             tetFaceI = fI;
 531
 532             label eIndex = -1;
 533
 534             if (edDir == 1)
 535             {
 536                 // Edge is in the forward circulation of this face, so
 537                 // work with the start point of the edge
 538
 539                 eIndex = findIndex(otherFace, e.start());
 540             }
 541             else
 542             {
 543                 // edDir == -1, so the edge is in the reverse
 544                 // circulation of this face, so work with the end
 545                 // point of the edge
 546
 547                 eIndex = findIndex(otherFace, e.end());
 548             }
 549
 550             label tetBasePtI = mesh_.tetBasePtIs()[fI];
 551
 552             if (tetBasePtI == -1)
 553             {
 554                 FatalErrorIn
 555                 (
 556                     "inline void "
 557                     "Foam::particle::crossEdgeConnectedFace"
 558                     "("
 559                         "const label& cellI,"
 560                         "label& tetFaceI,"
 561                         "label& tetPtI,"
 562                         "const edge& e"
 563                     ")"
 564                 )
 565                     << "No base point for face " << fI << ", " << f
 566                     << ", produces a decomposition that has a minimum "
 567                     << "volume greater than tolerance."
 568                     << abort(FatalError);
 569             }
 570
 571             // Find eIndex relative to the base point on new face
 572             eIndex -= tetBasePtI;
 573
 574             if (neg(eIndex))
 575             {
 576                 eIndex = (eIndex + otherFace.size()) % otherFace.size();
 577             }
 578
 579             if (eIndex == 0)
 580             {
 581                 // The point is the base point, so this is first tet
 582                 // in the face circulation
 583
 584                 tetPtI = 1;
 585             }
 586             else if (eIndex == otherFace.size() - 1)
 587             {
 588                 // The point is the last before the base point, so
 589                 // this is the last tet in the face circulation
 590
 591                 tetPtI = otherFace.size() - 2;
 592             }
 593             else
 594             {
 595                 tetPtI = eIndex;
 596             }
 597
 598             break;
 599         }
 600     }
 601 }
 602
 603
 604 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 605
 606 inline Foam::label Foam::particle::getNewParticleID() const
 607 {
 608     label id = particleCount_++;
 609
 610     if (id == labelMax)
 611     {
 612         WarningIn("particle::getNewParticleID() const")
 613             << "Particle counter has overflowed. This might cause problems"
 614             << " when reconstructing particle tracks." << endl;
 615     }
 616     return id;
 617 }
 618
 619
 620 inline const Foam::polyMesh& Foam::particle::mesh() const
 621 {
 622     return mesh_;
 623 }
 624
 625
 626 inline const Foam::vector& Foam::particle::position() const
 627 {
 628     return position_;
 629 }
 630
 631
 632 inline Foam::vector& Foam::particle::position()
 633 {
 634     return position_;
 635 }
 636
 637
 638 inline Foam::label Foam::particle::cell() const
 639 {
 640     return cellI_;
 641 }
 642
 643
 644 inline Foam::label& Foam::particle::cell()
 645 {
 646     return cellI_;
 647 }
 648
 649
 650 inline Foam::label Foam::particle::tetFace() const
 651 {
 652     return tetFaceI_;
 653 }
 654
 655
 656 inline Foam::label& Foam::particle::tetFace()
 657 {
 658     return tetFaceI_;
 659 }
 660
 661
 662 inline Foam::label Foam::particle::tetPt() const
 663 {
 664     return tetPtI_;
 665 }
 666
 667
 668 inline Foam::label& Foam::particle::tetPt()
 669 {
 670     return tetPtI_;
 671 }
 672
 673
 674 inline Foam::tetIndices Foam::particle::currentTetIndices() const
 675 {
 676     return tetIndices(cellI_, tetFaceI_, tetPtI_, mesh_);
 677 }
 678
 679
 680 inline Foam::tetPointRef Foam::particle::currentTet() const
 681 {
 682     return currentTetIndices().tet(mesh_);
 683 }
 684
 685
 686 inline Foam::vector Foam::particle::normal() const
 687 {
 688     return currentTetIndices().faceTri(mesh_).normal();
 689 }
 690
 691
 692 inline Foam::vector Foam::particle::oldNormal() const
 693 {
 694     return currentTetIndices().oldFaceTri(mesh_).normal();
 695 }
 696
 697
 698 inline Foam::label Foam::particle::face() const
 699 {
 700     return faceI_;
 701 }
 702
 703
 704 inline Foam::label& Foam::particle::face()
 705 {
 706     return faceI_;
 707 }
 708
 709
 710 inline void Foam::particle::initCellFacePt()
 711 {
 712     if (cellI_ == -1)
 713     {
 714         mesh_.findCellFacePt
 715         (
 716             position_,
 717             cellI_,
 718             tetFaceI_,
 719             tetPtI_
 720         );
 721
 722         if (cellI_ == -1)
 723         {
 724             FatalErrorIn("void Foam::particle::initCellFacePt()")
 725                 << "cell, tetFace and tetPt search failure at position "
 726                 << position_ << abort(FatalError);
 727         }
 728     }
 729     else
 730     {
 731         mesh_.findTetFacePt(cellI_, position_, tetFaceI_, tetPtI_);
 732
 733         if (tetFaceI_ == -1 || tetPtI_ == -1)
 734         {
 735             label oldCellI = cellI_;
 736
 737             mesh_.findCellFacePt
 738             (
 739                 position_,
 740                 cellI_,
 741                 tetFaceI_,
 742                 tetPtI_
 743             );
 744
 745             if (cellI_ == -1 || tetFaceI_ == -1 || tetPtI_ == -1)
 746             {
 747                 // The particle has entered this function with a cell
 748                 // number, but hasn't been able to find a cell to
 749                 // occupy.
 750
 751                 if(!mesh_.pointInCellBB(position_, oldCellI, 0.1))
 752                 {
 753                     // If the position is not inside the (slightly
 754                     // extended) bound-box of the cell that it thought
 755                     // it should be in, then this is considered an
 756                     // error.
 757
 758                     FatalErrorIn("void Foam::particle::initCellFacePt()")
 759                         << "cell, tetFace and tetPt search failure at position "
 760                         << position_ << nl << "for requested cell " << oldCellI
 761                         << abort(FatalError);
 762                 }
 763
 764                 // The position is in the (slightly extended)
 765                 // bound-box of the cell.  This situation may arise
 766                 // because the face decomposition of the cell is not
 767                 // the same as when the particle acquired the cell
 768                 // index.  For example, it has been read into a mesh
 769                 // that has made a different face base-point decision
 770                 // for a boundary face and now this particle is in a
 771                 // position that is not in the mesh.  Gradually move
 772                 // the particle towards the centre of the cell that it
 773                 // thought that it was in.
 774
 775                 cellI_ = oldCellI;
 776
 777                 point newPosition = position_;
 778
 779                 const point& cC = mesh_.cellCentres()[cellI_];
 780
 781                 label trap(1.0/trackingCorrectionTol + 1);
 782
 783                 label iterNo = 0;
 784
 785                 do
 786                 {
 787                     newPosition += trackingCorrectionTol*(cC - position_);
 788
 789                     mesh_.findTetFacePt
 790                     (
 791                         cellI_,
 792                         newPosition,
 793                         tetFaceI_,
 794                         tetPtI_
 795                     );
 796
 797                     iterNo++;
 798
 799                 } while (tetFaceI_ < 0  && iterNo <= trap);
 800
 801                 if (tetFaceI_ == -1)
 802                 {
 803                     FatalErrorIn("void Foam::particle::initCellFacePt()")
 804                         << "cell, tetFace and tetPt search failure at position "
 805                         << position_ << abort(FatalError);
 806                 }
 807
 808                 if (debug)
 809                 {
 810                     WarningIn("void Foam::particle::initCellFacePt()")
 811                         << "Particle moved from " << position_
 812                         << " to " << newPosition
 813                         << " in cell " << cellI_
 814                         << " tetFace " << tetFaceI_
 815                         << " tetPt " << tetPtI_ << nl
 816                         << "    (A fraction of "
 817                         << 1.0 - mag(cC - newPosition)/mag(cC - position_)
 818                         << " of the distance to the cell centre)"
 819                         << " because a decomposition tetFace and tetPt "
 820                         << "could not be found."
 821                         << endl;
 822                 }
 823
 824                 position_ = newPosition;
 825             }
 826
 827             if (debug && cellI_ != oldCellI)
 828             {
 829                 WarningIn("void Foam::particle::initCellFacePt()")
 830                     << "Particle at position " << position_
 831                     << " searched for a cell, tetFace and tetPt." << nl
 832                     << "    Found"
 833                     << " cell " << cellI_
 834                     << " tetFace " << tetFaceI_
 835                     << " tetPt " << tetPtI_ << nl
 836                     << "    This is a different cell to that which was supplied"
 837                     << " (" << oldCellI << ")." << nl
 838                     << endl;
 839             }
 840         }
 841     }
 842 }
 843
 844
 845 inline bool Foam::particle::onBoundary() const
 846 {
 847     return faceI_ != -1 && faceI_ >= mesh_.nInternalFaces();
 848 }
 849
 850
 851 inline Foam::scalar& Foam::particle::stepFraction()
 852 {
 853     return stepFraction_;
 854 }
 855
 856
 857 inline Foam::scalar Foam::particle::stepFraction() const
 858 {
 859     return stepFraction_;
 860 }
 861
 862
 863 inline Foam::label Foam::particle::origProc() const
 864 {
 865     return origProc_;
 866 }
 867
 868
 869 inline Foam::label& Foam::particle::origProc()
 870 {
 871     return origProc_;
 872 }
 873
 874
 875 inline Foam::label Foam::particle::origId() const
 876 {
 877     return origId_;
 878 }
 879
 880
 881 inline Foam::label& Foam::particle::origId()
 882 {
 883     return origId_;
 884 }
 885
 886
 887 inline bool Foam::particle::softImpact() const
 888 {
 889     return false;
 890 }
 891
 892
 893 inline Foam::scalar Foam::particle::currentTime() const
 894 {
 895     return
 896         mesh_.time().value()
 897       + stepFraction_*mesh_.time().deltaTValue();
 898 }
 899
 900
 901 inline bool Foam::particle::internalFace(const label faceI) const
 902 {
 903     return mesh_.isInternalFace(faceI);
 904 }
 905
 906
 907 bool Foam::particle::boundaryFace(const label faceI) const
 908 {
 909     return !internalFace(faceI);
 910 }
 911
 912
 913 inline Foam::label Foam::particle::patch(const label faceI) const
 914 {
 915     return mesh_.boundaryMesh().whichPatch(faceI);
 916 }
 917
 918
 919 inline Foam::label Foam::particle::patchFace
 920 (
 921     const label patchI,
 922     const label faceI
 923 ) const
 924 {
 925     return mesh_.boundaryMesh()[patchI].whichFace(faceI);
 926 }
 927
 928
 929 inline Foam::label Foam::particle::faceInterpolation() const
 930 {
 931     return faceI_;
 932 }
 933
 934
 935 // ************************************************************************* //