src/lagrangian/basic/particle/particleTemplates.C

   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 "IOPosition.H"
  27
  28 #include "cyclicPolyPatch.H"
  29 #include "processorPolyPatch.H"
  30 #include "symmetryPolyPatch.H"
  31 #include "wallPolyPatch.H"
  32 #include "wedgePolyPatch.H"
  33 #include "meshTools.H"
  34
  35 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  36
  37 template<class TrackData>
  38 void Foam::particle::prepareForParallelTransfer
  39 (
  40     const label patchI,
  41     TrackData& td
  42 )
  43 {
  44     // Convert the face index to be local to the processor patch
  45     faceI_ = patchFace(patchI, faceI_);
  46 }
  47
  48
  49 template<class TrackData>
  50 void Foam::particle::correctAfterParallelTransfer
  51 (
  52     const label patchI,
  53     TrackData& td
  54 )
  55 {
  56     const coupledPolyPatch& ppp =
  57         refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchI]);
  58
  59     cellI_ = ppp.faceCells()[faceI_];
  60
  61     // Have patch transform the position
  62     ppp.transformPosition(position_, faceI_);
  63
  64     // Transform the properties
  65     if (!ppp.parallel())
  66     {
  67         const tensor& T =
  68         (
  69             ppp.forwardT().size() == 1
  70           ? ppp.forwardT()[0]
  71           : ppp.forwardT()[faceI_]
  72         );
  73         transformProperties(T);
  74     }
  75     else if (ppp.separated())
  76     {
  77         const vector& s =
  78         (
  79             (ppp.separation().size() == 1)
  80           ? ppp.separation()[0]
  81           : ppp.separation()[faceI_]
  82         );
  83         transformProperties(-s);
  84     }
  85
  86     tetFaceI_ = faceI_ + ppp.start();
  87
  88     // Faces either side of a coupled patch have matched base indices,
  89     // tetPtI is specified relative to the base point, already and
  90     // opposite circulation directions by design, so if the vertices
  91     // are:
  92     // source:
  93     // face    (a b c d e f)
  94     // fPtI     0 1 2 3 4 5
  95     //            +
  96     // destination:
  97     // face    (a f e d c b)
  98     // fPtI     0 1 2 3 4 5
  99     //                  +
 100     // where a is the base point of the face are matching , and we
 101     // have fPtI = 1 on the source processor face, i.e. vertex b, then
 102     // this because of the face circulation direction change, vertex c
 103     // is the characterising point on the destination processor face,
 104     // giving the destination fPtI as:
 105     //     fPtI_d = f.size() - 1 - fPtI_s = 6 - 1 - 1 = 4
 106     // This relationship can be verified for other points and sizes of
 107     // face.
 108
 109     tetPtI_ = mesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
 110
 111     // Reset the face index for the next tracking operation
 112     if (stepFraction_ > (1.0 - SMALL))
 113     {
 114         stepFraction_ = 1.0;
 115         faceI_ = -1;
 116     }
 117     else
 118     {
 119         faceI_ += ppp.start();
 120     }
 121 }
 122
 123
 124 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 125
 126 template<class CloudType>
 127 void Foam::particle::readFields(CloudType& c)
 128 {
 129     if (!c.size())
 130     {
 131         return;
 132     }
 133
 134     IOobject procIO(c.fieldIOobject("origProcId", IOobject::MUST_READ));
 135
 136     if (procIO.headerOk())
 137     {
 138         IOField<label> origProcId(procIO);
 139         c.checkFieldIOobject(c, origProcId);
 140         IOField<label> origId(c.fieldIOobject("origId", IOobject::MUST_READ));
 141         c.checkFieldIOobject(c, origId);
 142
 143         label i = 0;
 144         forAllIter(typename CloudType, c, iter)
 145         {
 146             particle& p = iter();
 147
 148             p.origProc_ = origProcId[i];
 149             p.origId_ = origId[i];
 150             i++;
 151         }
 152     }
 153 }
 154
 155
 156 template<class CloudType>
 157 void Foam::particle::writeFields(const CloudType& c)
 158 {
 159     // Write the cloud position file
 160     IOPosition<CloudType> ioP(c);
 161     ioP.write();
 162
 163     label np =  c.size();
 164
 165     IOField<label> origProc
 166     (
 167         c.fieldIOobject("origProcId", IOobject::NO_READ),
 168         np
 169     );
 170     IOField<label> origId(c.fieldIOobject("origId", IOobject::NO_READ), np);
 171
 172     label i = 0;
 173     forAllConstIter(typename CloudType, c, iter)
 174     {
 175         origProc[i] = iter().origProc_;
 176         origId[i] = iter().origId_;
 177         i++;
 178     }
 179
 180     origProc.write();
 181     origId.write();
 182 }
 183
 184
 185 template<class TrackData>
 186 Foam::label Foam::particle::track(const vector& endPosition, TrackData& td)
 187 {
 188     faceI_ = -1;
 189
 190     // Tracks to endPosition or stop on boundary
 191     while (!onBoundary() && stepFraction_ < 1.0 - SMALL)
 192     {
 193         stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);
 194     }
 195
 196     return faceI_;
 197 }
 198
 199
 200 template<class TrackData>
 201 Foam::scalar Foam::particle::trackToFace
 202 (
 203     const vector& endPosition,
 204     TrackData& td
 205 )
 206 {
 207     typedef typename TrackData::cloudType cloudType;
 208     typedef typename cloudType::particleType particleType;
 209
 210     cloudType& cloud = td.cloud();
 211
 212
 213     const faceList& pFaces = mesh_.faces();
 214     const pointField& pPts = mesh_.points();
 215     const vectorField& pC = mesh_.cellCentres();
 216
 217     faceI_ = -1;
 218
 219     // Pout<< "Particle " << origId_ << " " << origProc_
 220     //     << " Tracking from " << position_
 221     //     << " to " << endPosition
 222     //     << endl;
 223
 224     // Pout<< "stepFraction " << stepFraction_ << nl
 225     //     << "cellI " << cellI_ << nl
 226     //     << "tetFaceI " << tetFaceI_ << nl
 227     //     << "tetPtI " << tetPtI_
 228     //     << endl;
 229
 230     scalar trackFraction = 0.0;
 231
 232     // Minimum tetrahedron decomposition of each cell of the mesh into
 233     // using the cell centre, base point on face, and further two
 234     // points on the face.  For each face of n points, there are n - 2
 235     // tets generated.
 236
 237     // The points for each tet are organised to match those used in the
 238     // tetrahedron class, supplying them in the order:
 239     //     Cc, basePt, pA, pB
 240     // where:
 241     //   + Cc is the cell centre;
 242     //   + basePt is the base point on the face;
 243     //   + pA and pB are the remaining points on the face, such that
 244     //     the circulation, {basePt, pA, pB} produces a positive
 245     //     normal by the right-hand rule.  pA and pB are chosen from
 246     //     tetPtI_ do accomplish this depending if the cell owns the
 247     //     face, tetPtI_ is the vertex that characterises the tet, and
 248     //     is the first vertex on the tet when circulating around the
 249     //     face. Therefore, the same tetPtI represents the same face
 250     //     triangle for both the owner and neighbour cell.
 251     //
 252     // Each tet has its four triangles represented in the same order:
 253     // 0) tri joining a tet to the tet across the face in next cell.
 254     //    This is the triangle opposite Cc.
 255     // 1) tri joining a tet to the tet that is in the same cell, but
 256     //    belongs to the face that shares the edge of the current face
 257     //    that doesn't contain basePt.  This is the triangle opposite
 258     //    basePt.
 259
 260     // 2) tri joining a tet to the tet that is in the same cell, but
 261     //    belongs to the face that shares the tet-edge (basePt - pB).
 262     //    This may be on the same face, or a different one.  This is
 263     //    the triangle opposite basePt.  This is the triangle opposite
 264     //    pA.
 265
 266     // 4) tri joining a tet to the tet that is in the same cell, but
 267     //    belongs to the face that shares the tet-edge (basePt - pA).
 268     //    This may be on the same face, or a different one.  This is
 269     //    the triangle opposite basePt.  This is the triangle opposite
 270     //    pA.
 271
 272     // Which tri (0..3) of the tet has been crossed
 273     label triI = -1;
 274
 275     // Determine which face was actually crossed.  lambdaMin < SMALL
 276     // is considered a trigger for a tracking correction towards the
 277     // current tet centre.
 278     scalar lambdaMin = VGREAT;
 279
 280     DynamicList<label>& tris = cloud.labels();
 281
 282     // Tet indices that will be set by hitWallFaces if a wall face is
 283     // to be hit, or are set when any wall tri of a tet is hit.
 284     // Carries the description of the tet on which the cell face has
 285     // been hit.  For the case of being set in hitWallFaces, this may
 286     // be a different tet to the one that the particle occupies.
 287     tetIndices faceHitTetIs;
 288
 289     // What tolerance is appropriate the minimum lambda numerator and
 290     // denominator for tracking in this cell.
 291     scalar lambdaDistanceTolerance =
 292         lambdaDistanceToleranceCoeff*mesh_.cellVolumes()[cellI_];
 293
 294     do
 295     {
 296         if (triI != -1)
 297         {
 298             // Change tet ownership because a tri face has been crossed
 299             tetNeighbour(triI);
 300         }
 301
 302         const Foam::face& f = pFaces[tetFaceI_];
 303
 304         bool own = (mesh_.faceOwner()[tetFaceI_] == cellI_);
 305
 306         label tetBasePtI = mesh_.tetBasePtIs()[tetFaceI_];
 307
 308         label basePtI = f[tetBasePtI];
 309
 310         label facePtI = (tetPtI_ + tetBasePtI) % f.size();
 311         label otherFacePtI = f.fcIndex(facePtI);
 312
 313         label fPtAI = -1;
 314         label fPtBI = -1;
 315
 316         if (own)
 317         {
 318             fPtAI = facePtI;
 319             fPtBI = otherFacePtI;
 320         }
 321         else
 322         {
 323             fPtAI = otherFacePtI;
 324             fPtBI = facePtI;
 325         }
 326
 327         tetPointRef tet
 328         (
 329             pC[cellI_],
 330             pPts[basePtI],
 331             pPts[f[fPtAI]],
 332             pPts[f[fPtBI]]
 333         );
 334
 335         if (lambdaMin < SMALL)
 336         {
 337             // Apply tracking correction towards tet centre
 338
 339             if (debug)
 340             {
 341                 Pout<< "tracking rescue using tetCentre from " << position();
 342             }
 343
 344             position_ += trackingCorrectionTol*(tet.centre() - position_);
 345
 346             if (debug)
 347             {
 348                 Pout<< " to " << position() << " due to "
 349                     << (tet.centre() - position_) << endl;
 350             }
 351
 352             cloud.trackingRescue();
 353
 354             return trackFraction;
 355         }
 356
 357         if (triI != -1 && mesh_.moving())
 358         {
 359             // Mesh motion requires stepFraction to be correct for
 360             // each tracking portion, so trackToFace must return after
 361             // every lambda calculation.
 362             return trackFraction;
 363         }
 364
 365         FixedList<vector, 4> tetAreas;
 366
 367         tetAreas[0] = tet.Sa();
 368         tetAreas[1] = tet.Sb();
 369         tetAreas[2] = tet.Sc();
 370         tetAreas[3] = tet.Sd();
 371
 372         FixedList<label, 4> tetPlaneBasePtIs;
 373
 374         tetPlaneBasePtIs[0] = basePtI;
 375         tetPlaneBasePtIs[1] = f[fPtAI];
 376         tetPlaneBasePtIs[2] = basePtI;
 377         tetPlaneBasePtIs[3] = basePtI;
 378
 379         findTris
 380         (
 381             endPosition,
 382             tris,
 383             tet,
 384             tetAreas,
 385             tetPlaneBasePtIs,
 386             lambdaDistanceTolerance
 387         );
 388
 389         // Reset variables for new track
 390         triI = -1;
 391         lambdaMin = VGREAT;
 392
 393         // Pout<< "tris " << tris << endl;
 394
 395         // Sets a value for lambdaMin and faceI_ if a wall face is hit
 396         // by the track.
 397         hitWallFaces
 398         (
 399             cloud,
 400             position_,
 401             endPosition,
 402             lambdaMin,
 403             faceHitTetIs
 404         );
 405
 406         // Did not hit any tet tri faces, and no wall face has been
 407         // found to hit.
 408         if (tris.empty() && faceI_ < 0)
 409         {
 410             position_ = endPosition;
 411
 412             return 1.0;
 413         }
 414         else
 415         {
 416             // Loop over all found tris and see if any of them find a
 417             // lambda value smaller than that found for a wall face.
 418             forAll(tris, i)
 419             {
 420                 label tI = tris[i];
 421
 422                 scalar lam = tetLambda
 423                 (
 424                     position_,
 425                     endPosition,
 426                     triI,
 427                     tetAreas[tI],
 428                     tetPlaneBasePtIs[tI],
 429                     cellI_,
 430                     tetFaceI_,
 431                     tetPtI_,
 432                     lambdaDistanceTolerance
 433                 );
 434
 435                 if (lam < lambdaMin)
 436                 {
 437                     lambdaMin = lam;
 438
 439                     triI = tI;
 440                 }
 441             }
 442         }
 443
 444         if (triI == 0)
 445         {
 446             // This must be a cell face crossing
 447             faceI_ = tetFaceI_;
 448
 449             // Set the faceHitTetIs to those for the current tet in case a
 450             // wall interaction is required with the cell face
 451             faceHitTetIs = tetIndices
 452             (
 453                 cellI_,
 454                 tetFaceI_,
 455                 tetBasePtI,
 456                 fPtAI,
 457                 fPtBI,
 458                 tetPtI_
 459             );
 460         }
 461         else if (triI > 0)
 462         {
 463             // A tri was found to be crossed before a wall face was hit (if any)
 464             faceI_ = -1;
 465         }
 466
 467         // Pout<< "track loop " << position_ << " " << endPosition << nl
 468         //     << "    " << cellI_
 469         //     << "    " << faceI_
 470         //     << " " << tetFaceI_
 471         //     << " " << tetPtI_
 472         //     << " " << triI
 473         //     << " " << lambdaMin
 474         //     << " " << trackFraction
 475         //     << endl;
 476
 477         // Pout<< "# Tracking loop tet "
 478         //     << origId_ << " " << origProc_<< nl
 479         //     << "# face: " << tetFaceI_ << nl
 480         //     << "# tetPtI: " << tetPtI_ << nl
 481         //     << "# tetBasePtI: " << mesh_.tetBasePtIs()[tetFaceI_] << nl
 482         //     << "# tet.mag(): " << tet.mag() << nl
 483         //     << "# tet.quality(): " << tet.quality()
 484         //     << endl;
 485
 486         // meshTools::writeOBJ(Pout, tet.a());
 487         // meshTools::writeOBJ(Pout, tet.b());
 488         // meshTools::writeOBJ(Pout, tet.c());
 489         // meshTools::writeOBJ(Pout, tet.d());
 490
 491         // Pout<< "f 1 3 2" << nl
 492         //     << "f 2 3 4" << nl
 493         //     << "f 1 4 3" << nl
 494         //     << "f 1 2 4" << endl;
 495
 496         // The particle can be 'outside' the tet.  This will yield a
 497         // lambda larger than 1, or smaller than 0.  For values < 0,
 498         // the particle travels away from the tet and we don't move
 499         // the particle, only change tet/cell.  For values larger than
 500         // 1, we move the particle to endPosition before the tet/cell
 501         // change.
 502         if (lambdaMin > SMALL)
 503         {
 504             if (lambdaMin <= 1.0)
 505             {
 506                 trackFraction += lambdaMin*(1 - trackFraction);
 507
 508                 position_ += lambdaMin*(endPosition - position_);
 509             }
 510             else
 511             {
 512                 trackFraction = 1.0;
 513
 514                 position_ = endPosition;
 515             }
 516         }
 517         else
 518         {
 519             // Set lambdaMin to zero to force a towards-tet-centre
 520             // correction.
 521             lambdaMin = 0.0;
 522         }
 523
 524     } while (faceI_ < 0);
 525
 526     particleType& p = static_cast<particleType&>(*this);
 527     p.hitFace(td);
 528
 529     if (internalFace(faceI_))
 530     {
 531         // Change tet ownership because a tri face has been crossed,
 532         // in general this is:
 533         //     tetNeighbour(triI);
 534         // but triI must be 0;
 535         // No modifications are required for triI = 0, no call required to
 536         //     tetNeighbour(0);
 537
 538         if (cellI_ == mesh_.faceOwner()[faceI_])
 539         {
 540             cellI_ = mesh_.faceNeighbour()[faceI_];
 541         }
 542         else if (cellI_ == mesh_.faceNeighbour()[faceI_])
 543         {
 544             cellI_ = mesh_.faceOwner()[faceI_];
 545         }
 546         else
 547         {
 548             FatalErrorIn("Particle::trackToFace(const vector&, TrackData&)")
 549                 << "addressing failure" << abort(FatalError);
 550         }
 551     }
 552     else
 553     {
 554         label origFaceI = faceI_;
 555         label patchI = patch(faceI_);
 556
 557         // No action taken for tetPtI_ for tetFaceI_ here, handled by
 558         // patch interaction call or later during processor transfer.
 559
 560         if
 561         (
 562             !p.hitPatch
 563             (
 564                 mesh_.boundaryMesh()[patchI],
 565                 td,
 566                 patchI,
 567                 trackFraction,
 568                 faceHitTetIs
 569             )
 570         )
 571         {
 572             // Did patch interaction model switch patches?
 573             if (faceI_ != origFaceI)
 574             {
 575                 patchI = patch(faceI_);
 576             }
 577
 578             const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 579
 580             if (isA<wedgePolyPatch>(patch))
 581             {
 582                 p.hitWedgePatch
 583                 (
 584                     static_cast<const wedgePolyPatch&>(patch), td
 585                 );
 586             }
 587             else if (isA<symmetryPolyPatch>(patch))
 588             {
 589                 p.hitSymmetryPatch
 590                 (
 591                     static_cast<const symmetryPolyPatch&>(patch), td
 592                 );
 593             }
 594             else if (isA<cyclicPolyPatch>(patch))
 595             {
 596                 p.hitCyclicPatch
 597                 (
 598                     static_cast<const cyclicPolyPatch&>(patch), td
 599                 );
 600             }
 601             else if (isA<processorPolyPatch>(patch))
 602             {
 603                 p.hitProcessorPatch
 604                 (
 605                     static_cast<const processorPolyPatch&>(patch), td
 606                 );
 607             }
 608             else if (isA<wallPolyPatch>(patch))
 609             {
 610                 p.hitWallPatch
 611                 (
 612                     static_cast<const wallPolyPatch&>(patch), td, faceHitTetIs
 613                 );
 614             }
 615             else
 616             {
 617                 p.hitPatch(patch, td);
 618             }
 619         }
 620     }
 621
 622     if (lambdaMin < SMALL)
 623     {
 624         // Apply tracking correction towards tet centre.
 625         // Generate current tet to find centre to apply correction.
 626
 627         tetPointRef tet = currentTet();
 628
 629         if (debug)
 630         {
 631             Pout<< "tracking rescue for lambdaMin:" << lambdaMin
 632                 << "from " << position();
 633         }
 634
 635         position_ += trackingCorrectionTol*(tet.centre() - position_);
 636
 637         if
 638         (
 639             cloud.hasWallImpactDistance()
 640          && !internalFace(faceHitTetIs.face())
 641          && cloud.cellHasWallFaces()[faceHitTetIs.cell()]
 642         )
 643         {
 644             const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 645
 646             label fI = faceHitTetIs.face();
 647
 648             label patchI = patches.patchID()[fI - mesh_.nInternalFaces()];
 649
 650             if (isA<wallPolyPatch>(patches[patchI]))
 651             {
 652                 // In the case of collision with a wall where there is
 653                 // a non-zero wallImpactDistance, it is possible for
 654                 // there to be a tracking correction required to bring
 655                 // the particle into the domain, but the position of
 656                 // the particle is further from the wall than the tet
 657                 // centre, in which case the normal correction can be
 658                 // counter-productive, i.e. pushes the particle
 659                 // further out of the domain.  In this case it is the
 660                 // position that hit the wall that is in need of a
 661                 // rescue correction.
 662
 663                 triPointRef wallTri = faceHitTetIs.faceTri(mesh_);
 664
 665                 tetPointRef wallTet = faceHitTetIs.tet(mesh_);
 666
 667                 vector nHat = wallTri.normal();
 668                 nHat /= mag(nHat);
 669
 670                 const scalar r = p.wallImpactDistance(nHat);
 671
 672                 // Removing (approximately) the wallTri normal
 673                 // component of the existing correction, to avoid the
 674                 // situation where the existing correction in the wall
 675                 // normal direction is larger towards the wall than
 676                 // the new correction is away from it.
 677                 position_ +=
 678                     trackingCorrectionTol
 679                    *(
 680                         (wallTet.centre() - (position_ + r*nHat))
 681                       - (nHat & (tet.centre() - position_))*nHat
 682                     );
 683             }
 684         }
 685
 686         if (debug)
 687         {
 688             Pout<< " to " << position() << endl;
 689         }
 690
 691         cloud.trackingRescue();
 692     }
 693
 694     return trackFraction;
 695 }
 696
 697
 698 template<class CloudType>
 699 void Foam::particle::hitWallFaces
 700 (
 701     const CloudType& cloud,
 702     const vector& from,
 703     const vector& to,
 704     scalar& lambdaMin,
 705     tetIndices& closestTetIs
 706 )
 707 {
 708     typedef typename CloudType::particleType particleType;
 709
 710     if (!(cloud.hasWallImpactDistance() && cloud.cellHasWallFaces()[cellI_]))
 711     {
 712         return;
 713     }
 714
 715     particleType& p = static_cast<particleType&>(*this);
 716
 717     const faceList& pFaces = mesh_.faces();
 718
 719     const Foam::cell& thisCell = mesh_.cells()[cellI_];
 720
 721     scalar lambdaDistanceTolerance =
 722         lambdaDistanceToleranceCoeff*mesh_.cellVolumes()[cellI_];
 723
 724     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 725
 726     forAll(thisCell, cFI)
 727     {
 728         label fI = thisCell[cFI];
 729
 730         if (internalFace(fI))
 731         {
 732             continue;
 733         }
 734
 735         label patchI = patches.patchID()[fI - mesh_.nInternalFaces()];
 736
 737         if (isA<wallPolyPatch>(patches[patchI]))
 738         {
 739             // Get the decomposition of this wall face
 740
 741             const List<tetIndices> faceTetIs =
 742                 polyMeshTetDecomposition::faceTetIndices(mesh_, fI, cellI_);
 743
 744             const Foam::face& f = pFaces[fI];
 745
 746             forAll(faceTetIs, tI)
 747             {
 748                 const tetIndices& tetIs = faceTetIs[tI];
 749
 750                 triPointRef tri = tetIs.faceTri(mesh_);
 751
 752                 vector n = tri.normal();
 753
 754                 vector nHat = n/mag(n);
 755
 756                 // Radius of particle with respect to this wall face
 757                 // triangle.  Assuming that the wallImpactDistance
 758                 // does not change as the particle or the mesh moves
 759                 // forward in time.
 760                 scalar r = p.wallImpactDistance(nHat);
 761
 762                 vector toPlusRNHat = to + r*nHat;
 763
 764                 // triI = 0 because it is the cell face tri of the tet
 765                 // we are concerned with.
 766                 scalar tetClambda = tetLambda
 767                 (
 768                     tetIs.tet(mesh_).centre(),
 769                     toPlusRNHat,
 770                     0,
 771                     n,
 772                     f[tetIs.faceBasePt()],
 773                     cellI_,
 774                     fI,
 775                     tetIs.tetPt(),
 776                     lambdaDistanceTolerance
 777                 );
 778
 779                 if ((tetClambda <= 0.0) || (tetClambda >= 1.0))
 780                 {
 781                     // toPlusRNHat is not on the outside of the plane of
 782                     // the wall face tri, the tri cannot be hit.
 783                     continue;
 784                 }
 785
 786                 // Check if the actual trajectory of the near-tri
 787                 // points intersects the triangle.
 788
 789                 vector fromPlusRNHat = from + r*nHat;
 790
 791                 // triI = 0 because it is the cell face tri of the tet
 792                 // we are concerned with.
 793                 scalar lambda = tetLambda
 794                 (
 795                     fromPlusRNHat,
 796                     toPlusRNHat,
 797                     0,
 798                     n,
 799                     f[tetIs.faceBasePt()],
 800                     cellI_,
 801                     fI,
 802                     tetIs.tetPt(),
 803                     lambdaDistanceTolerance
 804                 );
 805
 806                 pointHit hitInfo(vector::zero);
 807
 808                 if (mesh_.moving())
 809                 {
 810                     // For a moving mesh, the position of wall
 811                     // triangle needs to be moved in time to be
 812                     // consistent with the moment defined by the
 813                     // current value of stepFraction and the value of
 814                     // lambda just calculated.
 815
 816                     // Total fraction thought the timestep of the
 817                     // motion, including stepFraction before the
 818                     // current tracking step and the current
 819                     // lambda
 820                     // i.e.
 821                     // let s = stepFraction, l = lambda
 822                     // Motion of x in time:
 823                     // |-----------------|---------|---------|
 824                     // x00               x0        xi        x
 825                     //
 826                     // where xi is the correct value of x at the required
 827                     // tracking instant.
 828                     //
 829                     // x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
 830                     //
 831                     // i.e. the motion covered by previous tracking portions
 832                     // within this timestep, and
 833                     //
 834                     // xi = x0 + l*(x - x0)
 835                     //    = l*x + (1 - l)*x0
 836                     //    = l*x + (1 - l)*(s*x + (1 - s)*x00)
 837                     //    = (s + l - s*l)*x + (1 - (s + l - s*l))*x00
 838                     //
 839                     // let m = (s + l - s*l)
 840                     //
 841                     // xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
 842                     //
 843                     // In the same form as before.
 844
 845                     // Clip lambda to 0.0-1.0 to ensure that sensible
 846                     // positions are used for triangle intersections.
 847                     scalar lam = max(0.0, min(1.0, lambda));
 848
 849                     scalar m = stepFraction_ + lam - (stepFraction_*lam);
 850
 851                     triPointRef tri00 = tetIs.oldFaceTri(mesh_);
 852
 853                     // Use SMALL positive tolerance to make the triangle
 854                     // slightly "fat" to improve robustness.  Intersection
 855                     // is calculated as the ray (from + r*nHat) -> (to +
 856                     // r*nHat).
 857
 858                     point tPtA = tri00.a() + m*(tri.a() - tri00.a());
 859                     point tPtB = tri00.b() + m*(tri.b() - tri00.b());
 860                     point tPtC = tri00.c() + m*(tri.c() - tri00.c());
 861
 862                     triPointRef t(tPtA, tPtB, tPtC);
 863
 864                     // The point fromPlusRNHat + m*(to - from) is on the
 865                     // plane of the triangle.  Determine the
 866                     // intersection with this triangle by testing if
 867                     // this point is inside or outside of the triangle.
 868                     hitInfo = t.intersection
 869                     (
 870                         fromPlusRNHat + m*(to - from),
 871                         t.normal(),
 872                         intersection::FULL_RAY,
 873                         SMALL
 874                     );
 875                 }
 876                 else
 877                 {
 878                     // Use SMALL positive tolerance to make the triangle
 879                     // slightly "fat" to improve robustness.  Intersection
 880                     // is calculated as the ray (from + r*nHat) -> (to +
 881                     // r*nHat).
 882                     hitInfo = tri.intersection
 883                     (
 884                         fromPlusRNHat,
 885                         (to - from),
 886                         intersection::FULL_RAY,
 887                         SMALL
 888                     );
 889                 }
 890
 891                 if (hitInfo.hit())
 892                 {
 893                     if (lambda < lambdaMin)
 894                     {
 895                         lambdaMin = lambda;
 896
 897                         faceI_ = fI;
 898
 899                         closestTetIs = tetIs;
 900                     }
 901                 }
 902             }
 903         }
 904     }
 905 }
 906
 907
 908 template<class TrackData>
 909 void Foam::particle::hitFace(TrackData&)
 910 {}
 911
 912
 913 template<class TrackData>
 914 bool Foam::particle::hitPatch
 915 (
 916     const polyPatch&,
 917     TrackData&,
 918     const label,
 919     const scalar,
 920     const tetIndices&
 921 )
 922 {
 923     return false;
 924 }
 925
 926
 927 template<class TrackData>
 928 void Foam::particle::hitWedgePatch
 929 (
 930     const wedgePolyPatch& wpp,
 931     TrackData&
 932 )
 933 {
 934     FatalErrorIn
 935     (
 936         "void Foam::particle::hitWedgePatch"
 937         "("
 938             "const wedgePolyPatch& wpp, "
 939             "TrackData&"
 940         ")"
 941     )   << "Hitting a wedge patch should not be possible."
 942         << abort(FatalError);
 943
 944     vector nf = normal();
 945     nf /= mag(nf);
 946
 947     transformProperties(I - 2.0*nf*nf);
 948 }
 949
 950
 951 template<class TrackData>
 952 void Foam::particle::hitSymmetryPatch
 953 (
 954     const symmetryPolyPatch& spp,
 955     TrackData&
 956 )
 957 {
 958     vector nf = normal();
 959     nf /= mag(nf);
 960
 961     transformProperties(I - 2.0*nf*nf);
 962 }
 963
 964
 965 template<class TrackData>
 966 void Foam::particle::hitCyclicPatch
 967 (
 968     const cyclicPolyPatch& cpp,
 969     TrackData& td
 970 )
 971 {
 972     faceI_ = cpp.transformGlobalFace(faceI_);
 973
 974     cellI_ = mesh_.faceOwner()[faceI_];
 975
 976     tetFaceI_ = faceI_;
 977
 978     // See note in correctAfterParallelTransfer for tetPtI_ addressing.
 979     tetPtI_ = mesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
 980
 981     const cyclicPolyPatch& receiveCpp = cpp.neighbPatch();
 982     label patchFacei = receiveCpp.whichFace(faceI_);
 983
 984     // Now the particle is on the receiving side
 985
 986     // Have patch transform the position
 987     receiveCpp.transformPosition(position_, patchFacei);
 988
 989     // Transform the properties
 990     if (!receiveCpp.parallel())
 991     {
 992         const tensor& T =
 993         (
 994             receiveCpp.forwardT().size() == 1
 995           ? receiveCpp.forwardT()[0]
 996           : receiveCpp.forwardT()[patchFacei]
 997         );
 998         transformProperties(T);
 999     }
1000     else if (receiveCpp.separated())
1001     {
1002         const vector& s =
1003         (
1004             (receiveCpp.separation().size() == 1)
1005           ? receiveCpp.separation()[0]
1006           : receiveCpp.separation()[patchFacei]
1007         );
1008         transformProperties(-s);
1009     }
1010 }
1011
1012
1013 template<class TrackData>
1014 void Foam::particle::hitProcessorPatch(const processorPolyPatch&, TrackData&)
1015 {}
1016
1017
1018 template<class TrackData>
1019 void Foam::particle::hitWallPatch
1020 (
1021     const wallPolyPatch&,
1022     TrackData&,
1023     const tetIndices&
1024 )
1025 {}
1026
1027
1028 template<class TrackData>
1029 void Foam::particle::hitPatch(const polyPatch&, TrackData&)
1030 {}
1031
1032
1033 // ************************************************************************* //