particleI.H
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-2015 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  return 0.0;
118  }
119  else
120  {
121  if (mag((to - from)) < tol/mag(n))
122  {
123  // 'Zero' length track (compared to the tolerance, which is
124  // based on the cell volume, divided by the tet face area), not
125  // along the face, face cannot be crossed.
126  return GREAT;
127  }
128  else
129  {
130  // Trajectory is non-zero and parallel to face
131  lambdaDenominator = sign(lambdaDenominator)*SMALL;
132  }
133  }
134  }
135 
136  return lambdaNumerator/lambdaDenominator;
137 }
138 
139 
140 inline Foam::scalar Foam::particle::movingTetLambda
141 (
142  const vector& from,
143  const vector& to,
144  const label triI,
145  const vector& n,
146  const label tetPlaneBasePtI,
147  const label cellI,
148  const label tetFaceI,
149  const label tetPtI,
150  const scalar tol
151 ) const
152 {
153  const pointField& pPts = mesh_.points();
154  const pointField& oldPPts = mesh_.oldPoints();
155 
156  // Base point of plane at end of motion
157  const point& b = pPts[tetPlaneBasePtI];
158 
159  // n: Normal of plane at end of motion
160 
161  // Base point of plane at start of timestep
162  const point& b00 = oldPPts[tetPlaneBasePtI];
163 
164  // Base point of plane at start of tracking portion (cast forward by
165  // stepFraction)
166  point b0 = b00 + stepFraction_*(b - b00);
167 
168  // Normal of plane at start of tracking portion
169  vector n0 = vector::zero;
170 
171  {
172  tetIndices tetIs(cellI, tetFaceI, tetPtI, mesh_);
173 
174  // Tet at timestep start
175  tetPointRef tet00 = tetIs.oldTet(mesh_);
176 
177  // Tet at timestep end
178  tetPointRef tet = tetIs.tet(mesh_);
179 
180  point tet0PtA = tet00.a() + stepFraction_*(tet.a() - tet00.a());
181  point tet0PtB = tet00.b() + stepFraction_*(tet.b() - tet00.b());
182  point tet0PtC = tet00.c() + stepFraction_*(tet.c() - tet00.c());
183  point tet0PtD = tet00.d() + stepFraction_*(tet.d() - tet00.d());
184 
185  // Tracking portion start tet (cast forward by stepFraction)
186  tetPointRef tet0(tet0PtA, tet0PtB, tet0PtC, tet0PtD);
187 
188  switch (triI)
189  {
190  case 0:
191  {
192  n0 = tet0.Sa();
193  break;
194  }
195  case 1:
196  {
197  n0 = tet0.Sb();
198  break;
199  }
200  case 2:
201  {
202  n0 = tet0.Sc();
203  break;
204  }
205  case 3:
206  {
207  n0 = tet0.Sd();
208  break;
209  }
210  default:
211  {
212  break;
213  }
214  }
215  }
216 
217  if (mag(n0) < SMALL)
218  {
219  // If the old normal is zero (for example in layer addition)
220  // then use the current normal;
221  n0 = n;
222  }
223 
224  scalar lambdaNumerator = 0;
225  scalar lambdaDenominator = 0;
226 
227  vector dP = to - from;
228  vector dN = n - n0;
229  vector dB = b - b0;
230  vector dS = from - b0;
231 
232  if (mag(dN) > SMALL)
233  {
234  scalar a = (dP - dB) & dN;
235  scalar b = ((dP - dB) & n0) + (dS & dN);
236  scalar c = dS & n0;
237 
238  if (mag(a) > SMALL)
239  {
240 
241  // Solve quadratic for lambda
242  scalar discriminant = sqr(b) - 4.0*a*c;
243 
244  if (discriminant < 0)
245  {
246  // Imaginary roots only - face not crossed
247  return GREAT;
248  }
249  else
250  {
251  scalar q = -0.5*(b + sign(b)*Foam::sqrt(discriminant));
252 
253  if (mag(q) < VSMALL)
254  {
255  // If q is zero, then l1 = q/a is the required
256  // value of lambda, and is zero.
257  return 0.0;
258  }
259 
260  scalar l1 = q/a;
261  scalar l2 = c/q;
262 
263  // There will be two roots, a big one and a little
264  // one, choose the little one.
265 
266  if (mag(l1) < mag(l2))
267  {
268  return l1;
269  }
270  else
271  {
272  return l2;
273  }
274  }
275  }
276  {
277  // When a is zero, solve the first order polynomial
278  lambdaNumerator = -c;
279  lambdaDenominator = b;
280  }
281  }
282  else
283  {
284  // When n = n0 is zero, there is no plane rotation, solve the
285  // first order polynomial
286  lambdaNumerator = -(dS & n0);
287  lambdaDenominator = ((dP - dB) & n0);
288  }
289 
290  if (mag(lambdaDenominator) < tol)
291  {
292  if (mag(lambdaNumerator) < tol)
293  {
294  // Track starts on the face, and is potentially
295  // parallel to it. +-tol)/+-tol is not a good
296  // comparison, return 0.0, in anticipation of tet
297  // centre correction.
298  return 0.0;
299  }
300  else
301  {
302  if (mag((to - from)) < tol/mag(n))
303  {
304  // Zero length track, not along the face, face
305  // cannot be crossed.
306  return GREAT;
307  }
308  else
309  {
310  // Trajectory is non-zero and parallel to face
311  lambdaDenominator = sign(lambdaDenominator)*SMALL;
312  }
313  }
314  }
315 
316  return lambdaNumerator/lambdaDenominator;
317 }
318 
319 
320 
322 {
323  const labelList& pOwner = mesh_.faceOwner();
324  const faceList& pFaces = mesh_.faces();
325 
326  bool own = (pOwner[tetFaceI_] == cellI_);
327 
328  const Foam::face& f = pFaces[tetFaceI_];
329 
330  label tetBasePtI = mesh_.tetBasePtIs()[tetFaceI_];
331 
332  if (tetBasePtI == -1)
333  {
335  (
336  "inline void Foam::particle::tetNeighbour(label triI)"
337  )
338  << "No base point for face " << tetFaceI_ << ", " << f
339  << ", produces a valid tet decomposition."
340  << abort(FatalError);
341  }
342 
343  label facePtI = (tetPtI_ + tetBasePtI) % f.size();
344  label otherFacePtI = f.fcIndex(facePtI);
345 
346  switch (triI)
347  {
348  case 0:
349  {
350  // Crossing this triangle changes tet to that in the
351  // neighbour cell over tetFaceI
352 
353  // Modification of cellI_ will happen by other indexing,
354  // tetFaceI_ and tetPtI don't change.
355 
356  break;
357  }
358  case 1:
359  {
361  (
362  cellI_,
363  tetFaceI_,
364  tetPtI_,
365  Foam::edge(f[facePtI], f[otherFacePtI])
366  );
367 
368  break;
369  }
370  case 2:
371  {
372  if (own)
373  {
374  if (tetPtI_ < f.size() - 2)
375  {
376  tetPtI_ = f.fcIndex(tetPtI_);
377  }
378  else
379  {
381  (
382  cellI_,
383  tetFaceI_,
384  tetPtI_,
385  Foam::edge(f[tetBasePtI], f[otherFacePtI])
386  );
387  }
388  }
389  else
390  {
391  if (tetPtI_ > 1)
392  {
393  tetPtI_ = f.rcIndex(tetPtI_);
394  }
395  else
396  {
398  (
399  cellI_,
400  tetFaceI_,
401  tetPtI_,
402  Foam::edge(f[tetBasePtI], f[facePtI])
403  );
404  }
405  }
406 
407  break;
408  }
409  case 3:
410  {
411  if (own)
412  {
413  if (tetPtI_ > 1)
414  {
415  tetPtI_ = f.rcIndex(tetPtI_);
416  }
417  else
418  {
420  (
421  cellI_,
422  tetFaceI_,
423  tetPtI_,
424  Foam::edge(f[tetBasePtI], f[facePtI])
425  );
426  }
427  }
428  else
429  {
430  if (tetPtI_ < f.size() - 2)
431  {
432  tetPtI_ = f.fcIndex(tetPtI_);
433  }
434  else
435  {
437  (
438  cellI_,
439  tetFaceI_,
440  tetPtI_,
441  Foam::edge(f[tetBasePtI], f[otherFacePtI])
442  );
443  }
444  }
445 
446  break;
447  }
448  default:
449  {
451  (
452  "inline void "
453  "Foam::particle::tetNeighbour(label triI)"
454  )
455  << "Tet tri face index error, can only be 0..3, supplied "
456  << triI << abort(FatalError);
457 
458  break;
459  }
460  }
461 }
462 
463 
465 (
466  const label& cellI,
467  label& tetFaceI,
468  label& tetPtI,
469  const edge& e
470 )
471 {
472  const faceList& pFaces = mesh_.faces();
473  const cellList& pCells = mesh_.cells();
474 
475  const Foam::face& f = pFaces[tetFaceI];
476 
477  const Foam::cell& thisCell = pCells[cellI];
478 
479  forAll(thisCell, cFI)
480  {
481  // Loop over all other faces of this cell and
482  // find the one that shares this edge
483 
484  label fI = thisCell[cFI];
485 
486  if (tetFaceI == fI)
487  {
488  continue;
489  }
490 
491  const Foam::face& otherFace = pFaces[fI];
492 
493  label edDir = otherFace.edgeDirection(e);
494 
495  if (edDir == 0)
496  {
497  continue;
498  }
499  else if (f == pFaces[fI])
500  {
501  // This is a necessary condition if using duplicate baffles
502  // (so coincident faces). We need to make sure we don't cross into
503  // the face with the same vertices since we might enter a tracking
504  // loop where it never exits. This test should be cheap
505  // for most meshes so can be left in for 'normal' meshes.
506  continue;
507  }
508  else
509  {
510  //Found edge on other face
511  tetFaceI = fI;
512 
513  label eIndex = -1;
514 
515  if (edDir == 1)
516  {
517  // Edge is in the forward circulation of this face, so
518  // work with the start point of the edge
519  eIndex = findIndex(otherFace, e.start());
520  }
521  else
522  {
523  // edDir == -1, so the edge is in the reverse
524  // circulation of this face, so work with the end
525  // point of the edge
526  eIndex = findIndex(otherFace, e.end());
527  }
528 
529  label tetBasePtI = mesh_.tetBasePtIs()[fI];
530 
531  if (tetBasePtI == -1)
532  {
534  (
535  "inline void "
536  "Foam::particle::crossEdgeConnectedFace"
537  "("
538  "const label& cellI,"
539  "label& tetFaceI,"
540  "label& tetPtI,"
541  "const edge& e"
542  ")"
543  )
544  << "No base point for face " << fI << ", " << f
545  << ", produces a decomposition that has a minimum "
546  << "volume greater than tolerance."
547  << abort(FatalError);
548  }
549 
550  // Find eIndex relative to the base point on new face
551  eIndex -= tetBasePtI;
552 
553  if (neg(eIndex))
554  {
555  eIndex = (eIndex + otherFace.size()) % otherFace.size();
556  }
557 
558  if (eIndex == 0)
559  {
560  // The point is the base point, so this is first tet
561  // in the face circulation
562  tetPtI = 1;
563  }
564  else if (eIndex == otherFace.size() - 1)
565  {
566  // The point is the last before the base point, so
567  // this is the last tet in the face circulation
568  tetPtI = otherFace.size() - 2;
569  }
570  else
571  {
572  tetPtI = eIndex;
573  }
574 
575  break;
576  }
577  }
578 }
579 
580 
581 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
582 
584 {
585  label id = particleCount_++;
586 
587  if (id == labelMax)
588  {
589  WarningIn("particle::getNewParticleID() const")
590  << "Particle counter has overflowed. This might cause problems"
591  << " when reconstructing particle tracks." << endl;
592  }
593  return id;
594 }
595 
596 
598 {
599  return mesh_;
600 }
601 
602 
604 {
605  return position_;
606 }
607 
608 
610 {
611  return position_;
612 }
613 
614 
616 {
617  return cellI_;
618 }
619 
620 
622 {
623  return cellI_;
624 }
625 
626 
628 {
629  return tetFaceI_;
630 }
631 
632 
634 {
635  return tetFaceI_;
636 }
637 
638 
640 {
641  return tetPtI_;
642 }
643 
644 
646 {
647  return tetPtI_;
648 }
649 
650 
652 {
654 }
655 
656 
658 {
659  return currentTetIndices().tet(mesh_);
660 }
661 
662 
664 {
666 }
667 
668 
670 {
672 }
673 
674 
676 {
677  return faceI_;
678 }
679 
680 
682 {
683  return faceI_;
684 }
685 
686 
688 {
689  if (cellI_ == -1)
690  {
692  (
693  position_,
694  cellI_,
695  tetFaceI_,
696  tetPtI_
697  );
698 
699  if (cellI_ == -1)
700  {
701  FatalErrorIn("void Foam::particle::initCellFacePt()")
702  << "cell, tetFace and tetPt search failure at position "
703  << position_ << abort(FatalError);
704  }
705  }
706  else
707  {
709 
710  if (tetFaceI_ == -1 || tetPtI_ == -1)
711  {
712  label oldCellI = cellI_;
713 
715  (
716  position_,
717  cellI_,
718  tetFaceI_,
719  tetPtI_
720  );
721 
722  if (cellI_ == -1 || tetFaceI_ == -1 || tetPtI_ == -1)
723  {
724  // The particle has entered this function with a cell
725  // number, but hasn't been able to find a cell to
726  // occupy.
727 
728  if (!mesh_.pointInCellBB(position_, oldCellI, 0.1))
729  {
730  // If the position is not inside the (slightly
731  // extended) bound-box of the cell that it thought
732  // it should be in, then this is considered an
733  // error.
734 
735  FatalErrorIn("void Foam::particle::initCellFacePt()")
736  << " cell, tetFace and tetPt search failure at "
737  << "position " << position_ << nl
738  << " for requested cell " << oldCellI << nl
739  << " If this is a restart or "
740  "reconstruction/decomposition etc. it is likely that"
741  " the write precision is not sufficient.\n"
742  " Either increase 'writePrecision' or "
743  "set 'writeFormat' to 'binary'"
744  << abort(FatalError);
745  }
746 
747  // The position is in the (slightly extended)
748  // bound-box of the cell. This situation may arise
749  // because the face decomposition of the cell is not
750  // the same as when the particle acquired the cell
751  // index. For example, it has been read into a mesh
752  // that has made a different face base-point decision
753  // for a boundary face and now this particle is in a
754  // position that is not in the mesh. Gradually move
755  // the particle towards the centre of the cell that it
756  // thought that it was in.
757 
758  cellI_ = oldCellI;
759 
760  point newPosition = position_;
761 
762  const point& cC = mesh_.cellCentres()[cellI_];
763 
764  label trap(1.0/trackingCorrectionTol + 1);
765 
766  label iterNo = 0;
767 
768  do
769  {
770  newPosition += trackingCorrectionTol*(cC - position_);
771 
773  (
774  cellI_,
775  newPosition,
776  tetFaceI_,
777  tetPtI_
778  );
779 
780  iterNo++;
781 
782  } while (tetFaceI_ < 0 && iterNo <= trap);
783 
784  if (tetFaceI_ == -1)
785  {
786  FatalErrorIn("void Foam::particle::initCellFacePt()")
787  << "cell, tetFace and tetPt search failure at position "
788  << position_ << abort(FatalError);
789  }
790 
791  if (debug)
792  {
793  WarningIn("void Foam::particle::initCellFacePt()")
794  << "Particle moved from " << position_
795  << " to " << newPosition
796  << " in cell " << cellI_
797  << " tetFace " << tetFaceI_
798  << " tetPt " << tetPtI_ << nl
799  << " (A fraction of "
800  << 1.0 - mag(cC - newPosition)/mag(cC - position_)
801  << " of the distance to the cell centre)"
802  << " because a decomposition tetFace and tetPt "
803  << "could not be found."
804  << endl;
805  }
806 
807  position_ = newPosition;
808  }
809 
810  if (debug && cellI_ != oldCellI)
811  {
812  WarningIn("void Foam::particle::initCellFacePt()")
813  << "Particle at position " << position_
814  << " searched for a cell, tetFace and tetPt." << nl
815  << " Found"
816  << " cell " << cellI_
817  << " tetFace " << tetFaceI_
818  << " tetPt " << tetPtI_ << nl
819  << " This is a different cell to that which was supplied"
820  << " (" << oldCellI << ")." << nl
821  << endl;
822  }
823  }
824  }
825 }
826 
827 
828 inline bool Foam::particle::onBoundary() const
829 {
830  return faceI_ != -1 && faceI_ >= mesh_.nInternalFaces();
831 }
832 
833 
834 inline Foam::scalar& Foam::particle::stepFraction()
835 {
836  return stepFraction_;
837 }
838 
839 
840 inline Foam::scalar Foam::particle::stepFraction() const
841 {
842  return stepFraction_;
843 }
844 
845 
847 {
848  return origProc_;
849 }
850 
851 
853 {
854  return origProc_;
855 }
856 
857 
859 {
860  return origId_;
861 }
862 
863 
865 {
866  return origId_;
867 }
868 
869 
870 inline bool Foam::particle::softImpact() const
871 {
872  return false;
873 }
874 
875 
876 inline Foam::scalar Foam::particle::currentTime() const
877 {
878  return
879  mesh_.time().value()
881 }
882 
883 
884 inline bool Foam::particle::internalFace(const label faceI) const
885 {
886  return mesh_.isInternalFace(faceI);
887 }
888 
889 
890 bool Foam::particle::boundaryFace(const label faceI) const
891 {
892  return !internalFace(faceI);
893 }
894 
895 
896 inline Foam::label Foam::particle::patch(const label faceI) const
897 {
898  return mesh_.boundaryMesh().whichPatch(faceI);
899 }
900 
901 
903 (
904  const label patchI,
905  const label faceI
906 ) const
907 {
908  return mesh_.boundaryMesh()[patchI].whichFace(faceI);
909 }
910 
911 
913 {
914  return faceI_;
915 }
916 
917 
918 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
scalar currentTime() const
Return the particle current time.
Definition: particleI.H:876
label & cell()
Return current cell particle is in.
Definition: particleI.H:621
void findTetFacePt(const label celli, const point &p, label &tetFacei, label &tetPti) const
Find the tetFacei and tetPti for point p in celli.
Definition: polyMesh.C:1341
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
const Point & b() const
Definition: tetrahedronI.H:80
dimensionedScalar neg(const dimensionedScalar &ds)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
label & face()
Return current face particle is on otherwise -1.
Definition: particleI.H:681
dimensioned< scalar > mag(const dimensioned< Type > &)
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
labelList f(nPoints)
const Point & d() const
Definition: tetrahedronI.H:94
const labelList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:864
triPointRef oldFaceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the.
Definition: tetIndicesI.H:140
label getNewParticleID() const
Get unique particle creation id.
Definition: particleI.H:583
label origId() const
Return const access to the particle id on originating processor.
Definition: particleI.H:858
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:242
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:310
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
const cellList & cells() const
dimensionedScalar sign(const dimensionedScalar &ds)
static label particleCount_
Cumulative particle counter - used to provode unique ID.
Definition: particle.H:315
const vectorField & cellCentres() const
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:53
label patchFace(const label patchI, const label faceI) const
Which face of this patch is this particle on.
Definition: particleI.H:903
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
Point centre() const
Return centre (centroid)
Definition: tetrahedronI.H:164
label origId_
Local particle id on originating processor.
Definition: particle.H:164
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet from the.
Definition: tetIndicesI.H:66
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
label tetFaceI_
Index of the face that owns the decomposed tet that the.
Definition: particle.H:153
label n
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static const scalar trackingCorrectionTol
Fraction of distance to tet centre to move a particle to.
Definition: particle.H:322
#define WarningIn(functionName)
Report a warning using Foam::Warning.
scalar deltaTValue() const
Return time step value.
Definition: TimeState.H:100
tetPointRef oldTet(const polyMesh &mesh) const
Return the geometry corresponding to this tet from the.
Definition: tetIndicesI.H:84
vector normal() const
Return the normal of the tri on tetFaceI_ for the.
Definition: particleI.H:663
vector normal() const
Return vector normal.
Definition: triangleI.H:103
const vector & position() const
Return current particle position.
Definition: particleI.H:603
tetPointRef currentTet() const
Return the geometry of the current tet that the.
Definition: particleI.H:657
vector oldNormal() const
Return the normal of the tri on tetFaceI_ for the.
Definition: particleI.H:669
label cellI_
Index of the cell it is in.
Definition: particle.H:143
#define forAll(list, i)
Definition: UList.H:421
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
const Point & a() const
Return vertices.
Definition: tetrahedronI.H:73
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
label faceInterpolation() const
Return the index of the face to be used in the interpolation.
Definition: particleI.H:912
const Point & c() const
Definition: tetrahedronI.H:87
const polyMesh & mesh_
Reference to the polyMesh database.
Definition: particle.H:137
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
scalar stepFraction_
Fraction of time-step completed.
Definition: particle.H:149
A tetrahedron primitive.
Definition: tetrahedron.H:62
const polyMesh & mesh() const
Return the mesh database.
Definition: particleI.H:597
label & tetFace()
Return current tet face particle is in.
Definition: particleI.H:633
label origProc() const
Return const access to the originating processor id.
Definition: particleI.H:846
errorManip< error > abort(error &err)
Definition: errorManip.H:131
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1073
bool softImpact() const
Return the impact model to be used, soft or hard (default).
Definition: particleI.H:870
bool onBoundary() const
Is the particle on the boundary/(or outside the domain)?
Definition: particleI.H:828
label otherFace(const primitiveMesh &, const label cellI, const label faceI, const label edgeI)
Return face on cell using edgeI but not faceI. Throws error.
Definition: meshTools.C:554
label & tetPt()
Return current tet face particle is in.
Definition: particleI.H:645
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
scalar & stepFraction()
Return the fraction of time-step completed.
Definition: particleI.H:834
label nInternalFaces() const
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
int edgeDirection(const edge &) const
Return the edge direction on the face.
Definition: face.C:784
error FatalError
label patch(const label faceI) const
Which patch is particle on.
Definition: particleI.H:896
label end() const
Return end vertex label.
Definition: edgeI.H:92
bool pointInCellBB(const point &p, label celli, scalar inflationFraction=0) const
Return true if the point in the cell bounding box.
triPointRef faceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the.
Definition: tetIndicesI.H:109
bool boundaryFace(const label faceI) const
Is this global face a boundary face?
Definition: particleI.H:890
const dimensionedScalar c
Speed of light in a vacuum.
void crossEdgeConnectedFace(const label &cellI, label &tetFaceI, label &tetPtI, const edge &e)
Cross the from the given face across the given edge of the.
Definition: particleI.H:465
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
Definition: UListI.H:65
const Time & time() const
Return time.
label start() const
Return start vertex label.
Definition: edgeI.H:81
bool internalFace(const label faceI) const
Is this global face an internal face?
Definition: particleI.H:884
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1060
scalar movingTetLambda(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 a moving tri face.
Definition: particleI.H:141
label tetPtI_
Index of the point on the face that defines the decomposed.
Definition: particle.H:158
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
void initCellFacePt()
Check the stored cell value (setting if necessary) and.
Definition: particleI.H:687
dimensionedSymmTensor sqr(const dimensionedVector &dv)
label origProc_
Originating processor id.
Definition: particle.H:161
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:73
void tetNeighbour(label triI)
Modify the tet owner data by crossing triI.
Definition: particleI.H:321
tetIndices currentTetIndices() const
Return the indices of the current tet that the.
Definition: particleI.H:651
label faceI_
Face index if the particle is on a face otherwise -1.
Definition: particle.H:146
void findCellFacePt(const point &p, label &celli, label &tetFacei, label &tetPti) const
Find the cell, tetFacei and tetPti for point p.
Definition: polyMesh.C:1261
vector position_
Position of particle.
Definition: particle.H:140
const Type & value() const
Return const reference to value.
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
static const label labelMax
Definition: label.H:62