particleTemplates.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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 "cyclicAMIPolyPatch.H"
30 #include "processorPolyPatch.H"
31 #include "symmetryPlanePolyPatch.H"
32 #include "symmetryPolyPatch.H"
33 #include "wallPolyPatch.H"
34 #include "wedgePolyPatch.H"
35 #include "meshTools.H"
36 
37 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
38 
39 template<class TrackData>
41 (
42  const label patchi,
43  TrackData& td
44 )
45 {
46  // Convert the face index to be local to the processor patch
47  facei_ = patchFace(patchi, facei_);
48 }
49 
50 
51 template<class TrackData>
53 (
54  const label patchi,
55  TrackData& td
56 )
57 {
58  const coupledPolyPatch& ppp =
59  refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]);
60 
61  celli_ = ppp.faceCells()[facei_];
62 
63  // Have patch transform the position
64  ppp.transformPosition(position_, facei_);
65 
66  // Transform the properties
67  if (!ppp.parallel())
68  {
69  const tensor& T =
70  (
71  ppp.forwardT().size() == 1
72  ? ppp.forwardT()[0]
73  : ppp.forwardT()[facei_]
74  );
75  transformProperties(T);
76  }
77  else if (ppp.separated())
78  {
79  const vector& s =
80  (
81  (ppp.separation().size() == 1)
82  ? ppp.separation()[0]
83  : ppp.separation()[facei_]
84  );
85  transformProperties(-s);
86  }
87 
88  tetFacei_ = facei_ + ppp.start();
89 
90  // Faces either side of a coupled patch have matched base indices,
91  // tetPtI is specified relative to the base point, already and
92  // opposite circulation directions by design, so if the vertices
93  // are:
94  // source:
95  // face (a b c d e f)
96  // fPtI 0 1 2 3 4 5
97  // +
98  // destination:
99  // face (a f e d c b)
100  // fPtI 0 1 2 3 4 5
101  // +
102  // where a is the base point of the face are matching , and we
103  // have fPtI = 1 on the source processor face, i.e. vertex b, then
104  // this because of the face circulation direction change, vertex c
105  // is the characterising point on the destination processor face,
106  // giving the destination fPtI as:
107  // fPtI_d = f.size() - 1 - fPtI_s = 6 - 1 - 1 = 4
108  // This relationship can be verified for other points and sizes of
109  // face.
110 
111  tetPtI_ = mesh_.faces()[tetFacei_].size() - 1 - tetPtI_;
112 
113  // Reset the face index for the next tracking operation
114  if (stepFraction_ > (1.0 - SMALL))
115  {
116  stepFraction_ = 1.0;
117  facei_ = -1;
118  }
119  else
120  {
121  facei_ += ppp.start();
122  }
123 }
124 
125 
126 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
127 
128 template<class CloudType>
130 {
131  if (!c.size())
132  {
133  return;
134  }
135 
136  IOobject procIO(c.fieldIOobject("origProcId", IOobject::MUST_READ));
137 
138  if (procIO.headerOk())
139  {
140  IOField<label> origProcId(procIO);
141  c.checkFieldIOobject(c, origProcId);
144 
145  label i = 0;
146  forAllIter(typename CloudType, c, iter)
147  {
148  particle& p = iter();
149 
150  p.origProc_ = origProcId[i];
151  p.origId_ = origId[i];
152  i++;
153  }
154  }
155 }
156 
157 
158 template<class CloudType>
160 {
161  // Write the cloud position file
162  IOPosition<CloudType> ioP(c);
163  ioP.write();
164 
165  label np = c.size();
166 
168  (
169  c.fieldIOobject("origProcId", IOobject::NO_READ),
170  np
171  );
173 
174  label i = 0;
175  forAllConstIter(typename CloudType, c, iter)
176  {
177  origProc[i] = iter().origProc_;
178  origId[i] = iter().origId_;
179  i++;
180  }
181 
182  origProc.write();
183  origId.write();
184 }
185 
186 
187 template<class TrackData>
188 Foam::label Foam::particle::track(const vector& endPosition, TrackData& td)
189 {
190  facei_ = -1;
191 
192  // Tracks to endPosition or stop on boundary
193  while (!onBoundary() && stepFraction_ < 1.0 - SMALL)
194  {
195  stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);
196  }
197 
198  return facei_;
199 }
200 
201 
202 template<class TrackData>
203 Foam::scalar Foam::particle::trackToFace
204 (
205  const vector& endPosition,
206  TrackData& td
207 )
208 {
209  typedef typename TrackData::cloudType cloudType;
210  typedef typename cloudType::particleType particleType;
211 
212  cloudType& cloud = td.cloud();
213 
214 
215  const faceList& pFaces = mesh_.faces();
216  const pointField& pPts = mesh_.points();
217  const vectorField& pC = mesh_.cellCentres();
218 
219  facei_ = -1;
220 
221  // Pout<< "Particle " << origId_ << " " << origProc_
222  // << " Tracking from " << position_
223  // << " to " << endPosition
224  // << endl;
225 
226  // Pout<< "stepFraction " << stepFraction_ << nl
227  // << "celli " << celli_ << nl
228  // << "tetFacei " << tetFacei_ << nl
229  // << "tetPtI " << tetPtI_
230  // << endl;
231 
232  scalar trackFraction = 0.0;
233 
234  // Minimum tetrahedron decomposition of each cell of the mesh into
235  // using the cell centre, base point on face, and further two
236  // points on the face. For each face of n points, there are n - 2
237  // tets generated.
238 
239  // The points for each tet are organised to match those used in the
240  // tetrahedron class, supplying them in the order:
241  // Cc, basePt, pA, pB
242  // where:
243  // + Cc is the cell centre;
244  // + basePt is the base point on the face;
245  // + pA and pB are the remaining points on the face, such that
246  // the circulation, {basePt, pA, pB} produces a positive
247  // normal by the right-hand rule. pA and pB are chosen from
248  // tetPtI_ do accomplish this depending if the cell owns the
249  // face, tetPtI_ is the vertex that characterises the tet, and
250  // is the first vertex on the tet when circulating around the
251  // face. Therefore, the same tetPtI represents the same face
252  // triangle for both the owner and neighbour cell.
253  //
254  // Each tet has its four triangles represented in the same order:
255  // 0) tri joining a tet to the tet across the face in next cell.
256  // This is the triangle opposite Cc.
257  // 1) tri joining a tet to the tet that is in the same cell, but
258  // belongs to the face that shares the edge of the current face
259  // that doesn't contain basePt. This is the triangle opposite
260  // basePt.
261 
262  // 2) tri joining a tet to the tet that is in the same cell, but
263  // belongs to the face that shares the tet-edge (basePt - pB).
264  // This may be on the same face, or a different one. This is
265  // the triangle opposite basePt. This is the triangle opposite
266  // pA.
267 
268  // 4) tri joining a tet to the tet that is in the same cell, but
269  // belongs to the face that shares the tet-edge (basePt - pA).
270  // This may be on the same face, or a different one. This is
271  // the triangle opposite basePt. This is the triangle opposite
272  // pA.
273 
274  // Which tri (0..3) of the tet has been crossed
275  label triI = -1;
276 
277  // Determine which face was actually crossed. lambdaMin < SMALL
278  // is considered a trigger for a tracking correction towards the
279  // current tet centre.
280  scalar lambdaMin = VGREAT;
281 
282  DynamicList<label>& tris = cloud.labels();
283 
284  // Tet indices that will be set by hitWallFaces if a wall face is
285  // to be hit, or are set when any wall tri of a tet is hit.
286  // Carries the description of the tet on which the cell face has
287  // been hit. For the case of being set in hitWallFaces, this may
288  // be a different tet to the one that the particle occupies.
289  tetIndices faceHitTetIs;
290 
291  // What tolerance is appropriate the minimum lambda numerator and
292  // denominator for tracking in this cell.
293  scalar lambdaDistanceTolerance =
295 
296  do
297  {
298  if (triI != -1)
299  {
300  // Change tet ownership because a tri face has been crossed
301  tetNeighbour(triI);
302  }
303 
304  const Foam::face& f = pFaces[tetFacei_];
305 
306  bool own = (mesh_.faceOwner()[tetFacei_] == celli_);
307 
308  label tetBasePtI = mesh_.tetBasePtIs()[tetFacei_];
309 
310  label basePtI = f[tetBasePtI];
311 
312  label facePtI = (tetPtI_ + tetBasePtI) % f.size();
313  label otherFacePtI = f.fcIndex(facePtI);
314 
315  label fPtAI = -1;
316  label fPtBI = -1;
317 
318  if (own)
319  {
320  fPtAI = facePtI;
321  fPtBI = otherFacePtI;
322  }
323  else
324  {
325  fPtAI = otherFacePtI;
326  fPtBI = facePtI;
327  }
328 
329  tetPointRef tet
330  (
331  pC[celli_],
332  pPts[basePtI],
333  pPts[f[fPtAI]],
334  pPts[f[fPtBI]]
335  );
336 
337  if (lambdaMin < SMALL)
338  {
339  // Apply tracking correction towards tet centre
340 
341  if (debug)
342  {
343  Pout<< "tracking rescue using tetCentre from " << position();
344  }
345 
347 
348  if (debug)
349  {
350  Pout<< " to " << position() << " due to "
351  << (tet.centre() - position_) << endl;
352  }
353 
354  cloud.trackingRescue();
355 
356  return trackFraction;
357  }
358 
359  if (triI != -1 && mesh_.moving())
360  {
361  // Mesh motion requires stepFraction to be correct for
362  // each tracking portion, so trackToFace must return after
363  // every lambda calculation.
364  return trackFraction;
365  }
366 
367  FixedList<vector, 4> tetAreas;
368 
369  tetAreas[0] = tet.Sa();
370  tetAreas[1] = tet.Sb();
371  tetAreas[2] = tet.Sc();
372  tetAreas[3] = tet.Sd();
373 
374  FixedList<label, 4> tetPlaneBasePtIs;
375 
376  tetPlaneBasePtIs[0] = basePtI;
377  tetPlaneBasePtIs[1] = f[fPtAI];
378  tetPlaneBasePtIs[2] = basePtI;
379  tetPlaneBasePtIs[3] = basePtI;
380 
381  findTris
382  (
383  endPosition,
384  tris,
385  tet,
386  tetAreas,
387  tetPlaneBasePtIs,
388  lambdaDistanceTolerance
389  );
390 
391  // Reset variables for new track
392  triI = -1;
393  lambdaMin = VGREAT;
394 
395  // Pout<< "tris " << tris << endl;
396 
397  // Sets a value for lambdaMin and facei_ if a wall face is hit
398  // by the track.
400  (
401  cloud,
402  position_,
403  endPosition,
404  lambdaMin,
405  faceHitTetIs
406  );
407 
408  // Did not hit any tet tri faces, and no wall face has been
409  // found to hit.
410  if (tris.empty() && facei_ < 0)
411  {
412  position_ = endPosition;
413 
414  return 1.0;
415  }
416  else
417  {
418  // Loop over all found tris and see if any of them find a
419  // lambda value smaller than that found for a wall face.
420  forAll(tris, i)
421  {
422  label tI = tris[i];
423 
424  scalar lam = tetLambda
425  (
426  position_,
427  endPosition,
428  tI,
429  tetAreas[tI],
430  tetPlaneBasePtIs[tI],
431  celli_,
432  tetFacei_,
433  tetPtI_,
434  lambdaDistanceTolerance
435  );
436 
437  if (lam < lambdaMin)
438  {
439  lambdaMin = lam;
440 
441  triI = tI;
442  }
443  }
444  }
445 
446  if (triI == 0)
447  {
448  // This must be a cell face crossing
449  facei_ = tetFacei_;
450 
451  // Set the faceHitTetIs to those for the current tet in case a
452  // wall interaction is required with the cell face
453  faceHitTetIs = tetIndices
454  (
455  celli_,
456  tetFacei_,
457  tetBasePtI,
458  fPtAI,
459  fPtBI,
460  tetPtI_
461  );
462  }
463  else if (triI > 0)
464  {
465  // A tri was found to be crossed before a wall face was hit (if any)
466  facei_ = -1;
467  }
468 
469  // Pout<< "track loop " << position_ << " " << endPosition << nl
470  // << " " << celli_
471  // << " " << facei_
472  // << " " << tetFacei_
473  // << " " << tetPtI_
474  // << " " << triI
475  // << " " << lambdaMin
476  // << " " << trackFraction
477  // << endl;
478 
479  // Pout<< "# Tracking loop tet "
480  // << origId_ << " " << origProc_<< nl
481  // << "# face: " << tetFacei_ << nl
482  // << "# tetPtI: " << tetPtI_ << nl
483  // << "# tetBasePtI: " << mesh_.tetBasePtIs()[tetFacei_] << nl
484  // << "# tet.mag(): " << tet.mag() << nl
485  // << "# tet.quality(): " << tet.quality()
486  // << endl;
487 
488  // meshTools::writeOBJ(Pout, tet.a());
489  // meshTools::writeOBJ(Pout, tet.b());
490  // meshTools::writeOBJ(Pout, tet.c());
491  // meshTools::writeOBJ(Pout, tet.d());
492 
493  // Pout<< "f 1 3 2" << nl
494  // << "f 2 3 4" << nl
495  // << "f 1 4 3" << nl
496  // << "f 1 2 4" << endl;
497 
498  // The particle can be 'outside' the tet. This will yield a
499  // lambda larger than 1, or smaller than 0. For values < 0,
500  // the particle travels away from the tet and we don't move
501  // the particle, only change tet/cell. For values larger than
502  // 1, we move the particle to endPosition before the tet/cell
503  // change.
504  if (lambdaMin > SMALL)
505  {
506  if (lambdaMin <= 1.0)
507  {
508  trackFraction += lambdaMin*(1 - trackFraction);
509  position_ += lambdaMin*(endPosition - position_);
510  }
511  else
512  {
513  position_ = endPosition;
514 
515  return 1.0;
516  }
517  }
518  else
519  {
520  // Set lambdaMin to zero to force a towards-tet-centre
521  // correction.
522  lambdaMin = 0.0;
523  }
524 
525  } while (facei_ < 0);
526 
527  particleType& p = static_cast<particleType&>(*this);
528  p.hitFace(td);
529 
530  if (internalFace(facei_))
531  {
532  // Change tet ownership because a tri face has been crossed,
533  // in general this is:
534  // tetNeighbour(triI);
535  // but triI must be 0;
536  // No modifications are required for triI = 0, no call required to
537  // tetNeighbour(0);
538 
539  if (celli_ == mesh_.faceOwner()[facei_])
540  {
542  }
543  else if (celli_ == mesh_.faceNeighbour()[facei_])
544  {
546  }
547  else
548  {
550  << "addressing failure" << abort(FatalError);
551  }
552  }
553  else
554  {
555  label origFacei = facei_;
556  label patchi = patch(facei_);
557 
558  // No action taken for tetPtI_ for tetFacei_ here, handled by
559  // patch interaction call or later during processor transfer.
560 
561  if
562  (
563  !p.hitPatch
564  (
566  td,
567  patchi,
568  trackFraction,
569  faceHitTetIs
570  )
571  )
572  {
573  // Did patch interaction model switch patches?
574  if (facei_ != origFacei)
575  {
576  patchi = patch(facei_);
577  }
578 
580 
581  if (isA<wedgePolyPatch>(patch))
582  {
583  p.hitWedgePatch
584  (
585  static_cast<const wedgePolyPatch&>(patch), td
586  );
587  }
588  else if (isA<symmetryPlanePolyPatch>(patch))
589  {
590  p.hitSymmetryPlanePatch
591  (
592  static_cast<const symmetryPlanePolyPatch&>(patch), td
593  );
594  }
595  else if (isA<symmetryPolyPatch>(patch))
596  {
597  p.hitSymmetryPatch
598  (
599  static_cast<const symmetryPolyPatch&>(patch), td
600  );
601  }
602  else if (isA<cyclicPolyPatch>(patch))
603  {
604  p.hitCyclicPatch
605  (
606  static_cast<const cyclicPolyPatch&>(patch), td
607  );
608  }
609  else if (isA<cyclicAMIPolyPatch>(patch))
610  {
611  p.hitCyclicAMIPatch
612  (
613  static_cast<const cyclicAMIPolyPatch&>(patch),
614  td,
615  endPosition - position_
616  );
617  }
618  else if (isA<processorPolyPatch>(patch))
619  {
620  p.hitProcessorPatch
621  (
622  static_cast<const processorPolyPatch&>(patch), td
623  );
624  }
625  else if (isA<wallPolyPatch>(patch))
626  {
627  p.hitWallPatch
628  (
629  static_cast<const wallPolyPatch&>(patch), td, faceHitTetIs
630  );
631  }
632  else
633  {
634  p.hitPatch(patch, td);
635  }
636  }
637  }
638 
639  if (lambdaMin < SMALL)
640  {
641  // Apply tracking correction towards tet centre.
642  // Generate current tet to find centre to apply correction.
643 
644  tetPointRef tet = currentTet();
645 
646  if (debug)
647  {
648  Pout<< "tracking rescue for lambdaMin:" << lambdaMin
649  << "from " << position();
650  }
651 
653 
654  if
655  (
656  cloud.hasWallImpactDistance()
657  && !internalFace(faceHitTetIs.face())
658  && cloud.cellHasWallFaces()[faceHitTetIs.cell()]
659  )
660  {
662 
663  label fI = faceHitTetIs.face();
664 
665  label patchi = patches.patchID()[fI - mesh_.nInternalFaces()];
666 
667  if (isA<wallPolyPatch>(patches[patchi]))
668  {
669  // In the case of collision with a wall where there is
670  // a non-zero wallImpactDistance, it is possible for
671  // there to be a tracking correction required to bring
672  // the particle into the domain, but the position of
673  // the particle is further from the wall than the tet
674  // centre, in which case the normal correction can be
675  // counter-productive, i.e. pushes the particle
676  // further out of the domain. In this case it is the
677  // position that hit the wall that is in need of a
678  // rescue correction.
679 
680  triPointRef wallTri = faceHitTetIs.faceTri(mesh_);
681 
682  tetPointRef wallTet = faceHitTetIs.tet(mesh_);
683 
684  vector nHat = wallTri.normal();
685  nHat /= mag(nHat);
686 
687  const scalar r = p.wallImpactDistance(nHat);
688 
689  // Removing (approximately) the wallTri normal
690  // component of the existing correction, to avoid the
691  // situation where the existing correction in the wall
692  // normal direction is larger towards the wall than
693  // the new correction is away from it.
694  position_ +=
696  *(
697  (wallTet.centre() - (position_ + r*nHat))
698  - (nHat & (tet.centre() - position_))*nHat
699  );
700  }
701  }
702 
703  if (debug)
704  {
705  Pout<< " to " << position() << endl;
706  }
707 
708  cloud.trackingRescue();
709  }
710 
711  return trackFraction;
712 }
713 
714 
715 template<class CloudType>
717 (
718  const CloudType& cloud,
719  const vector& from,
720  const vector& to,
721  scalar& lambdaMin,
722  tetIndices& closestTetIs
723 )
724 {
725  typedef typename CloudType::particleType particleType;
726 
727  if (!(cloud.hasWallImpactDistance() && cloud.cellHasWallFaces()[celli_]))
728  {
729  return;
730  }
731 
732  particleType& p = static_cast<particleType&>(*this);
733 
734  const faceList& pFaces = mesh_.faces();
735 
736  const Foam::cell& thisCell = mesh_.cells()[celli_];
737 
738  scalar lambdaDistanceTolerance =
740 
742 
743  forAll(thisCell, cFI)
744  {
745  label fI = thisCell[cFI];
746 
747  if (internalFace(fI))
748  {
749  continue;
750  }
751 
752  label patchi = patches.patchID()[fI - mesh_.nInternalFaces()];
753 
754  if (isA<wallPolyPatch>(patches[patchi]))
755  {
756  // Get the decomposition of this wall face
757 
758  const List<tetIndices> faceTetIs =
760 
761  const Foam::face& f = pFaces[fI];
762 
763  forAll(faceTetIs, tI)
764  {
765  const tetIndices& tetIs = faceTetIs[tI];
766 
767  triPointRef tri = tetIs.faceTri(mesh_);
768 
769  vector n = tri.normal();
770 
771  vector nHat = n/mag(n);
772 
773  // Radius of particle with respect to this wall face
774  // triangle. Assuming that the wallImpactDistance
775  // does not change as the particle or the mesh moves
776  // forward in time.
777  scalar r = p.wallImpactDistance(nHat);
778 
779  vector toPlusRNHat = to + r*nHat;
780 
781  // triI = 0 because it is the cell face tri of the tet
782  // we are concerned with.
783  scalar tetClambda = tetLambda
784  (
785  tetIs.tet(mesh_).centre(),
786  toPlusRNHat,
787  0,
788  n,
789  f[tetIs.faceBasePt()],
790  celli_,
791  fI,
792  tetIs.tetPt(),
793  lambdaDistanceTolerance
794  );
795 
796  if ((tetClambda <= 0.0) || (tetClambda >= 1.0))
797  {
798  // toPlusRNHat is not on the outside of the plane of
799  // the wall face tri, the tri cannot be hit.
800  continue;
801  }
802 
803  // Check if the actual trajectory of the near-tri
804  // points intersects the triangle.
805 
806  vector fromPlusRNHat = from + r*nHat;
807 
808  // triI = 0 because it is the cell face tri of the tet
809  // we are concerned with.
810  scalar lambda = tetLambda
811  (
812  fromPlusRNHat,
813  toPlusRNHat,
814  0,
815  n,
816  f[tetIs.faceBasePt()],
817  celli_,
818  fI,
819  tetIs.tetPt(),
820  lambdaDistanceTolerance
821  );
822 
823  pointHit hitInfo(Zero);
824 
825  if (mesh_.moving())
826  {
827  // For a moving mesh, the position of wall
828  // triangle needs to be moved in time to be
829  // consistent with the moment defined by the
830  // current value of stepFraction and the value of
831  // lambda just calculated.
832 
833  // Total fraction thought the timestep of the
834  // motion, including stepFraction before the
835  // current tracking step and the current
836  // lambda
837  // i.e.
838  // let s = stepFraction, l = lambda
839  // Motion of x in time:
840  // |-----------------|---------|---------|
841  // x00 x0 xi x
842  //
843  // where xi is the correct value of x at the required
844  // tracking instant.
845  //
846  // x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
847  //
848  // i.e. the motion covered by previous tracking portions
849  // within this timestep, and
850  //
851  // xi = x0 + l*(x - x0)
852  // = l*x + (1 - l)*x0
853  // = l*x + (1 - l)*(s*x + (1 - s)*x00)
854  // = (s + l - s*l)*x + (1 - (s + l - s*l))*x00
855  //
856  // let m = (s + l - s*l)
857  //
858  // xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
859  //
860  // In the same form as before.
861 
862  // Clip lambda to 0.0-1.0 to ensure that sensible
863  // positions are used for triangle intersections.
864  scalar lam = max(0.0, min(1.0, lambda));
865 
866  scalar m = stepFraction_ + lam - (stepFraction_*lam);
867 
868  triPointRef tri00 = tetIs.oldFaceTri(mesh_);
869 
870  // Use SMALL positive tolerance to make the triangle
871  // slightly "fat" to improve robustness. Intersection
872  // is calculated as the ray (from + r*nHat) -> (to +
873  // r*nHat).
874 
875  point tPtA = tri00.a() + m*(tri.a() - tri00.a());
876  point tPtB = tri00.b() + m*(tri.b() - tri00.b());
877  point tPtC = tri00.c() + m*(tri.c() - tri00.c());
878 
879  triPointRef t(tPtA, tPtB, tPtC);
880 
881  // The point fromPlusRNHat + m*(to - from) is on the
882  // plane of the triangle. Determine the
883  // intersection with this triangle by testing if
884  // this point is inside or outside of the triangle.
885  hitInfo = t.intersection
886  (
887  fromPlusRNHat + m*(to - from),
888  t.normal(),
890  SMALL
891  );
892  }
893  else
894  {
895  // Use SMALL positive tolerance to make the triangle
896  // slightly "fat" to improve robustness. Intersection
897  // is calculated as the ray (from + r*nHat) -> (to +
898  // r*nHat).
899  hitInfo = tri.intersection
900  (
901  fromPlusRNHat,
902  (to - from),
904  SMALL
905  );
906  }
907 
908  if (hitInfo.hit())
909  {
910  if (lambda < lambdaMin)
911  {
912  lambdaMin = lambda;
913 
914  facei_ = fI;
915 
916  closestTetIs = tetIs;
917  }
918  }
919  }
920  }
921  }
922 }
923 
924 
925 template<class TrackData>
926 void Foam::particle::hitFace(TrackData&)
927 {}
928 
929 
930 template<class TrackData>
932 (
933  const polyPatch&,
934  TrackData&,
935  const label,
936  const scalar,
937  const tetIndices&
938 )
939 {
940  return false;
941 }
942 
943 
944 template<class TrackData>
946 (
947  const wedgePolyPatch& wpp,
948  TrackData&
949 )
950 {
952  << "Hitting a wedge patch should not be possible."
953  << abort(FatalError);
954 
955  vector nf = normal();
956  nf /= mag(nf);
957 
958  transformProperties(I - 2.0*nf*nf);
959 }
960 
961 
962 template<class TrackData>
964 (
965  const symmetryPlanePolyPatch& spp,
966  TrackData&
967 )
968 {
969  vector nf = normal();
970  nf /= mag(nf);
971 
972  transformProperties(I - 2.0*nf*nf);
973 }
974 
975 
976 template<class TrackData>
978 (
979  const symmetryPolyPatch& spp,
980  TrackData&
981 )
982 {
983  vector nf = normal();
984  nf /= mag(nf);
985 
986  transformProperties(I - 2.0*nf*nf);
987 }
988 
989 
990 template<class TrackData>
992 (
993  const cyclicPolyPatch& cpp,
994  TrackData& td
995 )
996 {
998 
1000 
1001  tetFacei_ = facei_;
1002 
1003  // See note in correctAfterParallelTransfer for tetPtI_ addressing.
1004  tetPtI_ = mesh_.faces()[tetFacei_].size() - 1 - tetPtI_;
1005 
1006  const cyclicPolyPatch& receiveCpp = cpp.neighbPatch();
1007  label patchFacei = receiveCpp.whichFace(facei_);
1008 
1009  // Now the particle is on the receiving side
1010 
1011  // Have patch transform the position
1012  receiveCpp.transformPosition(position_, patchFacei);
1013 
1014  // Transform the properties
1015  if (!receiveCpp.parallel())
1016  {
1017  const tensor& T =
1018  (
1019  receiveCpp.forwardT().size() == 1
1020  ? receiveCpp.forwardT()[0]
1021  : receiveCpp.forwardT()[patchFacei]
1022  );
1024  }
1025  else if (receiveCpp.separated())
1026  {
1027  const vector& s =
1028  (
1029  (receiveCpp.separation().size() == 1)
1030  ? receiveCpp.separation()[0]
1031  : receiveCpp.separation()[patchFacei]
1032  );
1033  transformProperties(-s);
1034  }
1035 }
1036 
1037 
1038 template<class TrackData>
1041  const cyclicAMIPolyPatch& cpp,
1042  TrackData& td,
1043  const vector& direction
1044 )
1045 {
1046  const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch();
1047 
1048  // Patch face index on sending side
1049  label patchFacei = facei_ - cpp.start();
1050 
1051  // Patch face index on receiving side - also updates position
1052  patchFacei = cpp.pointFace(patchFacei, direction, position_);
1053 
1054  if (patchFacei < 0)
1055  {
1057  << "Particle lost across " << cyclicAMIPolyPatch::typeName
1058  << " patches " << cpp.name() << " and " << receiveCpp.name()
1059  << " at position " << position_ << abort(FatalError);
1060  }
1061 
1062  // Convert face index into global numbering
1063  facei_ = patchFacei + receiveCpp.start();
1064 
1066 
1067  tetFacei_ = facei_;
1068 
1069  // See note in correctAfterParallelTransfer for tetPtI_ addressing.
1070  tetPtI_ = mesh_.faces()[tetFacei_].size() - 1 - tetPtI_;
1071 
1072  // Now the particle is on the receiving side
1073 
1074  // Transform the properties
1075  if (!receiveCpp.parallel())
1076  {
1077  const tensor& T =
1078  (
1079  receiveCpp.forwardT().size() == 1
1080  ? receiveCpp.forwardT()[0]
1081  : receiveCpp.forwardT()[patchFacei]
1082  );
1084  }
1085  else if (receiveCpp.separated())
1086  {
1087  const vector& s =
1088  (
1089  (receiveCpp.separation().size() == 1)
1090  ? receiveCpp.separation()[0]
1091  : receiveCpp.separation()[patchFacei]
1092  );
1093  transformProperties(-s);
1094  }
1095 }
1096 
1097 
1098 template<class TrackData>
1100 {}
1101 
1102 
1103 template<class TrackData>
1106  const wallPolyPatch&,
1107  TrackData&,
1108  const tetIndices&
1109 )
1110 {}
1111 
1112 
1113 template<class TrackData>
1114 void Foam::particle::hitPatch(const polyPatch&, TrackData&)
1115 {}
1116 
1117 
1118 // ************************************************************************* //
Symmetry patch for non-planar or multi-plane patches.
ParcelType particleType
Definition: Cloud.H:114
void hitSymmetryPlanePatch(const symmetryPlanePolyPatch &, TrackData &td)
Overridable function to handle the particle hitting a.
bool onBoundary() const
Is the particle on the boundary/(or outside the domain)?
Definition: particleI.H:810
const labelList & patchID() const
Per boundary face label the patch index.
const cyclicPolyPatch & neighbPatch() const
vector Sa() const
Return face normal.
Definition: tetrahedronI.H:136
A tetrahedron primitive.
Definition: tetrahedron.H:62
void findTris(const vector &position, DynamicList< label > &faceList, const tetPointRef &tet, const FixedList< vector, 4 > &tetAreas, const FixedList< label, 4 > &tetPlaneBasePtIs, const scalar tol) const
Find the tet tri faces between position and tet centre.
Definition: particleI.H:32
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:57
label cell() const
Return the cell.
Definition: tetIndicesI.H:30
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
static const scalar trackingCorrectionTol
Fraction of distance to tet centre to move a particle to.
Definition: particle.H:326
virtual void transformPosition(pointField &l) const
Transform a patch-based position from other side to this side.
label whichFace(const label l) const
Return label of face in patch from global face label.
Definition: polyPatch.H:377
void prepareForParallelTransfer(const label patchi, TrackData &td)
Convert global addressing to the processor patch.
void checkFieldIOobject(const Cloud< ParticleType > &c, const IOField< DataType > &data) const
Check lagrangian data field.
Definition: CloudIO.C:215
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:53
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
virtual const vectorField & separation() const
If the planes are separated the separation vector.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:313
const Point & c() const
Return third vertex.
Definition: triangleI.H:82
triPointRef faceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the.
Definition: tetIndicesI.H:109
scalar tetLambda(const vector &from, const vector &to, const label triI, const vector &tetArea, const label tetPlaneBasePtI, const label celli, const label tetFacei, const label tetPtI, const scalar tol) const
Find the lambda value for the line to-from across the.
Definition: particleI.H:69
void hitWallPatch(const wallPolyPatch &, TrackData &td, const tetIndices &tetIs)
Overridable function to handle the particle hitting a wallPatch.
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:453
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
label tetPtI_
Index of the point on the face that defines the decomposed.
Definition: particle.H:158
void hitWallFaces(const CloudType &td, const vector &from, const vector &to, scalar &lambdaMin, tetIndices &closestTetIs)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
label track(const vector &endPosition, TrackData &td)
Track particle to end of trajectory.
static const scalar lambdaDistanceToleranceCoeff
Fraction of the cell volume to use in determining tolerance values.
Definition: particle.H:330
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
patches[0]
static void readFields(CloudType &c)
Read the fields associated with the owner cloud.
virtual void transformProperties(const tensor &T)
Transform the physical properties of the particle.
Definition: particle.C:123
vector position_
Position of particle.
Definition: particle.H:140
bool hitPatch(const polyPatch &, TrackData &td, const label patchi, const scalar trackFraction, const tetIndices &tetIs)
Overridable function to handle the particle hitting a.
void hitCyclicAMIPatch(const cyclicAMIPolyPatch &, TrackData &td, const vector &direction)
Overridable function to handle the particle hitting a cyclicAMIPatch.
label tetPt() const
Return the characterising tetPtI.
Definition: tetIndicesI.H:60
label origId() const
Return const access to the particle id on originating processor.
Definition: particleI.H:840
const cellList & cells() const
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
bool moving() const
Is mesh moving.
Definition: polyMesh.H:493
label origId_
Local particle id on originating processor.
Definition: particle.H:164
Base particle class.
Definition: particle.H:78
bool internalFace(const label facei) const
Is this global face an internal face?
Definition: particleI.H:866
label transformGlobalFace(const label facei) const
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
Neighbour processor patch.
label faceBasePt() const
Return the face base point.
Definition: tetIndicesI.H:42
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
tetPointRef currentTet() const
Return the geometry of the current tet that the.
Definition: particleI.H:640
label origProc_
Originating processor id.
Definition: particle.H:161
const vector & position() const
Return current particle position.
Definition: particleI.H:586
void hitProcessorPatch(const processorPolyPatch &, TrackData &td)
Overridable function to handle the particle hitting a.
virtual bool separated() const
Are the planes separated.
label patch(const label facei) const
Which patch is particle on.
Definition: particleI.H:878
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
Point centre() const
Return centre (centroid)
Definition: tetrahedronI.H:164
static void writeFields(const CloudType &c)
Write the fields associated with the owner cloud.
static const Identity< scalar > I
Definition: Identity.H:93
Cyclic plane patch.
Foam::wallPolyPatch.
Definition: wallPolyPatch.H:48
const Point & a() const
Return first vertex.
Definition: triangleI.H:70
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints,-1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){pointMap[i]=i;}for(label i=0;i< nPoints;i++){if(f[i] > 0.0){hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei){if(edges[ei].mag(points)< SMALL){label start=pointMap[edges[ei].start()];while(start!=pointMap[start]){start=pointMap[start];}label end=pointMap[edges[ei].end()];while(end!=pointMap[end]){end=pointMap[end];}label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;}}cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){cellShape &cs=cellShapes[celli];forAll(cs, i){cs[i]=pointMap[cs[i]];}cs.collapse();}label bcIDs[11]={-1, 0, 2, 4,-1, 5,-1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&symmetryPolyPatch::typeName,&wedgePolyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&symmetryPolyPatch::typeName,&oldCyclicPolyPatch::typeName};enum patchTypeNames{PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={"piston","valve","liner","cylinderHead","axis","wedge","inflow","outflow","presin","presout","symmetryPlane","cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
const polyMesh & mesh_
Reference to the polyMesh database.
Definition: particle.H:137
Wedge front and back plane patch.
triPointRef oldFaceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the.
Definition: tetIndicesI.H:140
const vectorField & cellCentres() const
static const zero Zero
Definition: zero.H:91
vector normal() const
Return the normal of the tri on tetFacei_ for the.
Definition: particleI.H:646
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
virtual bool parallel() const
Are the cyclic planes parallel.
Cyclic patch for Arbitrary Mesh Interface (AMI)
label face() const
Return the face.
Definition: tetIndicesI.H:36
void hitWedgePatch(const wedgePolyPatch &, TrackData &td)
Overridable function to handle the particle hitting a wedgePatch.
const word & name() const
Return name.
scalar stepFraction_
Fraction of time-step completed.
Definition: particle.H:149
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void hitFace(TrackData &td)
Overridable function to handle the particle hitting a face.
label celli_
Index of the cell it is in.
Definition: particle.H:143
vector Sb() const
Definition: tetrahedronI.H:143
Foam::polyBoundaryMesh.
label tetFacei_
Index of the face that owns the decomposed tet that the.
Definition: particle.H:153
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
label origProc() const
Return const access to the originating processor id.
Definition: particleI.H:828
const labelList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:818
labelList f(nPoints)
label size() const
Definition: Cloud.H:175
void tetNeighbour(label triI)
Modify the tet owner data by crossing triI.
Definition: particleI.H:321
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const scalarField & cellVolumes() const
const volScalarField & T
scalar trackToFace(const vector &endPosition, TrackData &td)
Track particle to a given position and returns 1.0 if the.
const Point & b() const
Return second vertex.
Definition: triangleI.H:76
label patchi
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1023
virtual bool write() const
Write using setting from DB.
Definition: IOPosition.C:51
virtual void transformPosition(pointField &) const =0
Transform a patch-based position from other side to this side.
const labelUList & faceCells() const
Return face-cell addressing.
Definition: polyPatch.C:315
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
label pointFace(const label facei, const vector &n, point &p) const
Return face index on neighbour patch which shares point p.
const dimensionedScalar c
Speed of light in a vacuum.
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
virtual const tensorField & forwardT() const
Return face transformation tensor.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
IOobject fieldIOobject(const word &fieldName, const IOobject::readOption r) const
Helper to construct IOobject for field and current time.
Definition: CloudIO.C:195
void correctAfterParallelTransfer(const label patchi, TrackData &td)
Convert processor patch addressing to the global equivalents.
virtual bool hasWallImpactDistance() const
Switch to specify if particles of the cloud can return.
Definition: Cloud.H:209
pointHit intersection(const point &p, const vector &q, const intersection::algorithm alg, const scalar tol=0.0) const
Fast intersection with a ray.
Definition: triangleI.H:425
volScalarField & p
const PackedBoolList & cellHasWallFaces() const
Whether each cell has any wall faces (demand driven data)
Definition: Cloud.C:149
Helper IO class to read and write particle positions.
Definition: Cloud.H:55
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1017
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
vector Sd() const
Definition: tetrahedronI.H:157
vector Sc() const
Definition: tetrahedronI.H:150
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
label nInternalFaces() const
A primitive field of type <T> with automated input and output.
Definition: IOField.H:50
void hitCyclicPatch(const cyclicPolyPatch &, TrackData &td)
Overridable function to handle the particle hitting a cyclicPatch.
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet from the.
Definition: tetIndicesI.H:66
cloud(const objectRegistry &, const word &cloudName="")
Construct for the given objectRegistry and named cloud instance.
Definition: cloud.C:42
virtual const cyclicAMIPolyPatch & neighbPatch() const
Return a reference to the neighbour patch.
void hitSymmetryPatch(const symmetryPolyPatch &, TrackData &td)
Overridable function to handle the particle hitting a.
vector normal() const
Return vector normal.
Definition: triangleI.H:103
static List< tetIndices > faceTetIndices(const polyMesh &mesh, label fI, label cI)
Return the tet decomposition of the given face, with.
label facei_
Face index if the particle is on a face otherwise -1.
Definition: particle.H:146