particle.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 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 "particle.H"
27 #include "polyTopoChangeMap.H"
28 #include "transform.H"
29 #include "treeDataCell.H"
30 #include "cubicEqn.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 const Foam::label Foam::particle::maxNTracksBehind_ = 48;
35 
36 Foam::label Foam::particle::particleCount_ = 0;
37 
38 namespace Foam
39 {
40  defineTypeNameAndDebug(particle, 0);
41 }
42 
43 
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 
46 void Foam::particle::stationaryTetReverseTransform
47 (
48  vector& centre,
49  scalar& detA,
51 ) const
52 {
53  barycentricTensor A = stationaryTetTransform();
54 
55  const vector ab = A.b() - A.a();
56  const vector ac = A.c() - A.a();
57  const vector ad = A.d() - A.a();
58  const vector bc = A.c() - A.b();
59  const vector bd = A.d() - A.b();
60 
61  centre = A.a();
62 
63  detA = ab & (ac ^ ad);
64 
66  (
67  bd ^ bc,
68  ac ^ ad,
69  ad ^ ab,
70  ab ^ ac
71  );
72 }
73 
74 
75 void Foam::particle::movingTetReverseTransform
76 (
77  const scalar fraction,
78  Pair<vector>& centre,
79  FixedList<scalar, 4>& detA,
80  FixedList<barycentricTensor, 3>& T
81 ) const
82 {
83  Pair<barycentricTensor> A = movingTetTransform(fraction);
84 
85  const Pair<vector> ab(A[0].b() - A[0].a(), A[1].b() - A[1].a());
86  const Pair<vector> ac(A[0].c() - A[0].a(), A[1].c() - A[1].a());
87  const Pair<vector> ad(A[0].d() - A[0].a(), A[1].d() - A[1].a());
88  const Pair<vector> bc(A[0].c() - A[0].b(), A[1].c() - A[1].b());
89  const Pair<vector> bd(A[0].d() - A[0].b(), A[1].d() - A[1].b());
90 
91  centre[0] = A[0].a();
92  centre[1] = A[1].a();
93 
94  detA[0] = ab[0] & (ac[0] ^ ad[0]);
95  detA[1] =
96  (ab[1] & (ac[0] ^ ad[0]))
97  + (ab[0] & (ac[1] ^ ad[0]))
98  + (ab[0] & (ac[0] ^ ad[1]));
99  detA[2] =
100  (ab[0] & (ac[1] ^ ad[1]))
101  + (ab[1] & (ac[0] ^ ad[1]))
102  + (ab[1] & (ac[1] ^ ad[0]));
103  detA[3] = ab[1] & (ac[1] ^ ad[1]);
104 
105  T[0] = barycentricTensor
106  (
107  bd[0] ^ bc[0],
108  ac[0] ^ ad[0],
109  ad[0] ^ ab[0],
110  ab[0] ^ ac[0]
111  );
112  T[1] = barycentricTensor
113  (
114  (bd[0] ^ bc[1]) + (bd[1] ^ bc[0]),
115  (ac[0] ^ ad[1]) + (ac[1] ^ ad[0]),
116  (ad[0] ^ ab[1]) + (ad[1] ^ ab[0]),
117  (ab[0] ^ ac[1]) + (ab[1] ^ ac[0])
118  );
119  T[2] = barycentricTensor
120  (
121  bd[1] ^ bc[1],
122  ac[1] ^ ad[1],
123  ad[1] ^ ab[1],
124  ab[1] ^ ac[1]
125  );
126 }
127 
128 
129 void Foam::particle::reflect()
130 {
131  Swap(coordinates_.c(), coordinates_.d());
132 }
133 
134 
135 void Foam::particle::rotate(const bool reverse)
136 {
137  if (!reverse)
138  {
139  scalar temp = coordinates_.b();
140  coordinates_.b() = coordinates_.c();
141  coordinates_.c() = coordinates_.d();
142  coordinates_.d() = temp;
143  }
144  else
145  {
146  scalar temp = coordinates_.d();
147  coordinates_.d() = coordinates_.c();
148  coordinates_.c() = coordinates_.b();
149  coordinates_.b() = temp;
150  }
151 }
152 
153 
154 void Foam::particle::changeTet(const label tetTriI)
155 {
156  if (debug)
157  {
158  Info << "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
159  }
160 
161  const bool isOwner = mesh_.faceOwner()[tetFacei_] == celli_;
162 
163  const label firstTetPtI = 1;
164  const label lastTetPtI = mesh_.faces()[tetFacei_].size() - 2;
165 
166  if (tetTriI == 1)
167  {
168  changeFace(tetTriI);
169  }
170  else if (tetTriI == 2)
171  {
172  if (isOwner)
173  {
174  if (tetPti_ == lastTetPtI)
175  {
176  changeFace(tetTriI);
177  }
178  else
179  {
180  reflect();
181  tetPti_ += 1;
182  }
183  }
184  else
185  {
186  if (tetPti_ == firstTetPtI)
187  {
188  changeFace(tetTriI);
189  }
190  else
191  {
192  reflect();
193  tetPti_ -= 1;
194  }
195  }
196  }
197  else if (tetTriI == 3)
198  {
199  if (isOwner)
200  {
201  if (tetPti_ == firstTetPtI)
202  {
203  changeFace(tetTriI);
204  }
205  else
206  {
207  reflect();
208  tetPti_ -= 1;
209  }
210  }
211  else
212  {
213  if (tetPti_ == lastTetPtI)
214  {
215  changeFace(tetTriI);
216  }
217  else
218  {
219  reflect();
220  tetPti_ += 1;
221  }
222  }
223  }
224  else
225  {
227  << "Changing tet without changing cell should only happen when the "
228  << "track is on triangle 1, 2 or 3."
229  << exit(FatalError);
230  }
231 }
232 
233 
234 void Foam::particle::changeFace(const label tetTriI)
235 {
236  if (debug)
237  {
238  Info << "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
239  }
240 
241  // Get the old topology
242  const triFace triOldIs(currentTetIndices().faceTriIs(mesh_));
243 
244  // Get the shared edge and the pre-rotation
245  edge sharedEdge;
246  if (tetTriI == 1)
247  {
248  sharedEdge = edge(triOldIs[1], triOldIs[2]);
249  }
250  else if (tetTriI == 2)
251  {
252  sharedEdge = edge(triOldIs[2], triOldIs[0]);
253  }
254  else if (tetTriI == 3)
255  {
256  sharedEdge = edge(triOldIs[0], triOldIs[1]);
257  }
258  else
259  {
261  << "Changing face without changing cell should only happen when the"
262  << " track is on triangle 1, 2 or 3."
263  << exit(FatalError);
264 
265  sharedEdge = edge(-1, -1);
266  }
267 
268  // Find the face in the same cell that shares the edge, and the
269  // corresponding tetrahedra point
270  tetPti_ = -1;
271  forAll(mesh_.cells()[celli_], cellFaceI)
272  {
273  const label newFaceI = mesh_.cells()[celli_][cellFaceI];
274  const class face& newFace = mesh_.faces()[newFaceI];
275  const label newOwner = mesh_.faceOwner()[newFaceI];
276 
277  // Exclude the current face
278  if (tetFacei_ == newFaceI)
279  {
280  continue;
281  }
282 
283  // Loop over the edges, looking for the shared one
284  const label edgeComp = newOwner == celli_ ? -1 : +1;
285  label edgeI = 0;
286  for
287  (
288  ;
289  edgeI < newFace.size()
290  && edge::compare(sharedEdge, newFace.faceEdge(edgeI)) != edgeComp;
291  ++ edgeI
292  );
293 
294  // If the face does not contain the edge, then move on to the next face
295  if (edgeI >= newFace.size())
296  {
297  continue;
298  }
299 
300  // Make the edge index relative to the base point
301  const label newBaseI = max(0, mesh_.tetBasePtIs()[newFaceI]);
302  edgeI = (edgeI - newBaseI + newFace.size()) % newFace.size();
303 
304  // If the edge is next the base point (i.e., the index is 0 or n - 1),
305  // then we swap it for the adjacent edge. This new edge is opposite the
306  // base point, and defines the tet with the original edge in it.
307  edgeI = min(max(1, edgeI), newFace.size() - 2);
308 
309  // Set the new face and tet point
310  tetFacei_ = newFaceI;
311  tetPti_ = edgeI;
312 
313  // Exit the loop now that the tet point has been found
314  break;
315  }
316 
317  if (tetPti_ == -1)
318  {
320  << "The search for an edge-connected face and tet-point failed."
321  << exit(FatalError);
322  }
323 
324  // Pre-rotation puts the shared edge opposite the base of the tetrahedron
325  if (sharedEdge.otherVertex(triOldIs[1]) == -1)
326  {
327  rotate(false);
328  }
329  else if (sharedEdge.otherVertex(triOldIs[2]) == -1)
330  {
331  rotate(true);
332  }
333 
334  // Get the new topology
335  const triFace triNewIs = currentTetIndices().faceTriIs(mesh_);
336 
337  // Reflect to account for the change of triangle orientation on the new face
338  reflect();
339 
340  // Post rotation puts the shared edge back in the correct location
341  if (sharedEdge.otherVertex(triNewIs[1]) == -1)
342  {
343  rotate(true);
344  }
345  else if (sharedEdge.otherVertex(triNewIs[2]) == -1)
346  {
347  rotate(false);
348  }
349 }
350 
351 
352 void Foam::particle::changeCell()
353 {
354  if (debug)
355  {
356  Info << "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
357  }
358 
359  // Set the cell to be the one on the other side of the face
360  const label ownerCellI = mesh_.faceOwner()[tetFacei_];
361  const bool isOwner = celli_ == ownerCellI;
362  celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
363 
364  // Reflect to account for the change of triangle orientation in the new cell
365  reflect();
366 }
367 
368 
369 void Foam::particle::locate
370 (
371  const vector& position,
372  label celli,
373  const bool boundaryFail,
374  const string boundaryMsg
375 )
376 {
377  if (debug)
378  {
379  Info << "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
380  }
381 
382  // Find the cell, if it has not been given
383  if (celli < 0)
384  {
385  celli = mesh_.cellTree().findInside(position);
386  }
387  if (celli < 0)
388  {
390  << "Cell not found for particle position " << position << "."
391  << exit(FatalError);
392  }
393  celli_ = celli;
394 
395  // Track from the centre of the cell to the desired position
396  const vector displacement = position - mesh_.cellCentres()[celli_];
397 
398  // Loop all cell tets to find the one containing the position. Track
399  // through each tet from the cell centre. If a tet contains the position
400  // then the track will end with a single trackToTri.
401  const class cell& c = mesh_.cells()[celli_];
402  scalar minF = vGreat;
403  label minTetFacei = -1, minTetPti = -1;
404  forAll(c, cellTetFacei)
405  {
406  const class face& f = mesh_.faces()[c[cellTetFacei]];
407  for (label tetPti = 1; tetPti < f.size() - 1; ++ tetPti)
408  {
409  coordinates_ = barycentric(1, 0, 0, 0);
410  tetFacei_ = c[cellTetFacei];
411  tetPti_ = tetPti;
412  facei_ = -1;
413  reset(1);
414 
415  label tetTriI = -1;
416  const scalar f = trackToTri(displacement, 0, tetTriI);
417 
418  if (tetTriI == -1)
419  {
420  return;
421  }
422 
423  if (f < minF)
424  {
425  minF = f;
426  minTetFacei = tetFacei_;
427  minTetPti = tetPti_;
428  }
429  }
430  }
431 
432  // The particle must be (hopefully only slightly) outside the cell. Track
433  // into the tet which got the furthest.
434  coordinates_ = barycentric(1, 0, 0, 0);
435  tetFacei_ = minTetFacei;
436  tetPti_ = minTetPti;
437  facei_ = -1;
438  reset(1);
439 
440  track(displacement, 0);
441  if (!onFace())
442  {
443  return;
444  }
445 
446  // If we are here then we hit a boundary
447  if (boundaryFail)
448  {
449  FatalErrorInFunction << boundaryMsg << exit(FatalError);
450  }
451  else
452  {
453  static label nWarnings = 0;
454  static const label maxNWarnings = 100;
455  if (nWarnings < maxNWarnings)
456  {
457  WarningInFunction << boundaryMsg.c_str() << endl;
458  ++ nWarnings;
459  }
460  if (nWarnings == maxNWarnings)
461  {
463  << "Suppressing any further warnings about particles being "
464  << "located outside of the mesh." << endl;
465  ++ nWarnings;
466  }
467  }
468 }
469 
470 
471 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
472 
474 (
475  const polyMesh& mesh,
476  const barycentric& coordinates,
477  const label celli,
478  const label tetFacei,
479  const label tetPti
480 )
481 :
482  mesh_(mesh),
483  coordinates_(coordinates),
484  celli_(celli),
485  tetFacei_(tetFacei),
486  tetPti_(tetPti),
487  facei_(-1),
488  stepFraction_(1),
489  stepFractionBehind_(0),
490  nTracksBehind_(0),
491  origProc_(Pstream::myProcNo()),
492  origId_(getNewParticleID())
493 {}
494 
495 
497 (
498  const polyMesh& mesh,
499  const vector& position,
500  const label celli
501 )
502 :
503  mesh_(mesh),
504  coordinates_(- vGreat, - vGreat, - vGreat, - vGreat),
505  celli_(celli),
506  tetFacei_(-1),
507  tetPti_(-1),
508  facei_(-1),
509  stepFraction_(1),
510  stepFractionBehind_(0),
511  nTracksBehind_(0),
512  origProc_(Pstream::myProcNo()),
513  origId_(getNewParticleID())
514 {
515  locate
516  (
517  position,
518  celli,
519  false,
520  "Particle initialised with a location outside of the mesh."
521  );
522 }
523 
524 
526 :
527  mesh_(p.mesh_),
528  coordinates_(p.coordinates_),
529  celli_(p.celli_),
530  tetFacei_(p.tetFacei_),
531  tetPti_(p.tetPti_),
532  facei_(p.facei_),
533  stepFraction_(p.stepFraction_),
534  stepFractionBehind_(p.stepFractionBehind_),
535  nTracksBehind_(p.nTracksBehind_),
536  origProc_(p.origProc_),
537  origId_(p.origId_)
538 {}
539 
540 
542 :
543  mesh_(mesh),
544  coordinates_(p.coordinates_),
545  celli_(p.celli_),
546  tetFacei_(p.tetFacei_),
547  tetPti_(p.tetPti_),
548  facei_(p.facei_),
549  stepFraction_(p.stepFraction_),
550  stepFractionBehind_(p.stepFractionBehind_),
551  nTracksBehind_(p.nTracksBehind_),
552  origProc_(p.origProc_),
553  origId_(p.origId_)
554 {}
555 
556 
557 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
558 
559 Foam::scalar Foam::particle::track
560 (
561  const vector& displacement,
562  const scalar fraction
563 )
564 {
565  if (debug)
566  {
567  Info << "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
568  }
569 
570  scalar f = trackToFace(displacement, fraction);
571 
572  while (onInternalFace())
573  {
574  changeCell();
575 
576  f *= trackToFace(f*displacement, f*fraction);
577  }
578 
579  return f;
580 }
581 
582 
583 Foam::scalar Foam::particle::trackToCell
584 (
585  const vector& displacement,
586  const scalar fraction
587 )
588 {
589  if (debug)
590  {
591  Info << "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
592  }
593 
594  const scalar f = trackToFace(displacement, fraction);
595 
596  if (onInternalFace())
597  {
598  changeCell();
599  }
600 
601  return f;
602 }
603 
604 
605 Foam::scalar Foam::particle::trackToFace
606 (
607  const vector& displacement,
608  const scalar fraction
609 )
610 {
611  if (debug)
612  {
613  Info << "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
614  }
615 
616  scalar f = 1;
617 
618  label tetTriI = onFace() ? 0 : -1;
619 
620  facei_ = -1;
621 
622  // Loop the tets in the current cell
623  while (nTracksBehind_ < maxNTracksBehind_)
624  {
625  f *= trackToTri(f*displacement, f*fraction, tetTriI);
626 
627  if (tetTriI == -1)
628  {
629  // The track has completed within the current tet
630  return 0;
631  }
632  else if (tetTriI == 0)
633  {
634  // The track has hit a face, so set the current face and return
635  facei_ = tetFacei_;
636  return f;
637  }
638  else
639  {
640  // Move to the next tet and continue the track
641  changeTet(tetTriI);
642  }
643  }
644 
645  // Warn if stuck, and incorrectly advance the step fraction to completion
647  << "Particle #" << origId_ << " got stuck at " << position() << endl;
648 
649  stepFraction_ += f*fraction;
650 
651  stepFractionBehind_ = 0;
652  nTracksBehind_ = 0;
653 
654  return 0;
655 }
656 
657 
659 (
660  const vector& displacement,
661  const scalar fraction,
662  label& tetTriI
663 )
664 {
665  const vector x0 = position();
666  const vector x1 = displacement;
667  const barycentric y0 = coordinates_;
668 
669  if (debug)
670  {
671  Info<< "Particle " << origId() << endl << "Tracking from " << x0
672  << " along " << x1 << " to " << x0 + x1 << endl;
673  }
674 
675  // Get the tet geometry
676  vector centre;
677  scalar detA;
679  stationaryTetReverseTransform(centre, detA, T);
680 
681  if (debug)
682  {
683  vector o, b, v1, v2;
684  stationaryTetGeometry(o, b, v1, v2);
685  Info<< "Tet points o=" << o << ", b=" << b
686  << ", v1=" << v1 << ", v2=" << v2 << endl
687  << "Tet determinant = " << detA << endl
688  << "Start local coordinates = " << y0 << endl;
689  }
690 
691  // Calculate the local tracking displacement
692  barycentric Tx1(x1 & T);
693 
694  if (debug)
695  {
696  Info<< "Local displacement = " << Tx1 << "/" << detA << endl;
697  }
698 
699  // Calculate the hit fraction
700  label iH = -1;
701  scalar muH = detA > vSmall ? 1/detA : vGreat;
702  for (label i = 0; i < 4; ++ i)
703  {
704  if (Tx1[i] < - vSmall && Tx1[i] < - mag(detA)*small)
705  {
706  scalar mu = - y0[i]/Tx1[i];
707 
708  if (debug)
709  {
710  Info<< "Hit on tet face " << i << " at local coordinate "
711  << y0 + mu*Tx1 << ", " << mu*detA*100 << "% of the "
712  << "way along the track" << endl;
713  }
714 
715  if (0 <= mu && mu < muH)
716  {
717  iH = i;
718  muH = mu;
719  }
720  }
721  }
722 
723  // If there has been no hit on a degenerate or inverted tet then the
724  // displacement must be within the round off error. Advance the step
725  // fraction without moving and return.
726  if (iH == -1 && muH == vGreat)
727  {
728  stepFraction_ += fraction;
729  return 0;
730  }
731 
732  // Set the new coordinates
733  barycentric yH = y0 + muH*Tx1;
734 
735  // Clamp to zero any negative coordinates generated by round-off error
736  for (label i = 0; i < 4; ++ i)
737  {
738  yH.replace(i, i == iH ? 0 : max(0, yH[i]));
739  }
740 
741  // Re-normalise if within the tet
742  if (iH == -1)
743  {
744  yH /= cmptSum(yH);
745  }
746 
747  // Set the new position and hit index
748  coordinates_ = yH;
749  tetTriI = iH;
750 
751  // Set the proportion of the track that has been completed
752  stepFraction_ += fraction*muH*detA;
753 
754  if (debug)
755  {
756  if (iH != -1)
757  {
758  Info<< "Track hit tet face " << iH << " first" << endl;
759  }
760  else
761  {
762  Info<< "Track hit no tet faces" << endl;
763  }
764  Info<< "End local coordinates = " << yH << endl
765  << "End global coordinates = " << position() << endl
766  << "Tracking displacement = " << position() - x0 << endl
767  << muH*detA*100 << "% of the step from "
768  << stepFraction_ - fraction*muH*detA << " to "
769  << stepFraction_ - fraction*muH*detA + fraction
770  << " completed" << endl << endl;
771  }
772 
773  // Accumulate fraction behind
774  if (muH*detA < small || nTracksBehind_ > 0)
775  {
776  stepFractionBehind_ += (fraction != 0 ? fraction : 1)*muH*detA;
777 
778  if (stepFractionBehind_ > rootSmall)
779  {
780  stepFractionBehind_ = 0;
781  nTracksBehind_ = 0;
782  }
783  else
784  {
785  ++ nTracksBehind_;
786  }
787  }
788 
789  return iH != -1 ? 1 - muH*detA : 0;
790 }
791 
792 
794 (
795  const vector& displacement,
796  const scalar fraction,
797  label& tetTriI
798 )
799 {
800  const vector x0 = position();
801  const vector x1 = displacement;
802  const barycentric y0 = coordinates_;
803 
804  if (debug)
805  {
806  Info<< "Particle " << origId() << endl << "Tracking from " << x0
807  << " along " << x1 << " to " << x0 + x1 << endl;
808  }
809 
810  // Get the tet geometry
811  Pair<vector> centre;
814  movingTetReverseTransform(fraction, centre, detA, T);
815 
816  if (debug)
817  {
818  Pair<vector> o, b, v1, v2;
819  movingTetGeometry(fraction, o, b, v1, v2);
820  Info<< "Tet points o=" << o[0] << ", b=" << b[0]
821  << ", v1=" << v1[0] << ", v2=" << v2[0] << endl
822  << "Tet determinant = " << detA[0] << endl
823  << "Start local coordinates = " << y0[0] << endl;
824  }
825 
826  // Get the relative global position
827  const vector x0Rel = x0 - centre[0];
828  const vector x1Rel = x1 - centre[1];
829 
830  // Form the determinant and hit equations
831  cubicEqn detAEqn(sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
832  const barycentric yC(1, 0, 0, 0);
833  const barycentric hitEqnA =
834  ((x1Rel & T[2]) + detA[3]*yC)*sqr(detA[0]);
835  const barycentric hitEqnB =
836  ((x1Rel & T[1]) + (x0Rel & T[2]) + detA[2]*yC)*detA[0];
837  const barycentric hitEqnC =
838  ((x1Rel & T[0]) + (x0Rel & T[1]) + detA[1]*yC);
839  const barycentric hitEqnD = y0;
840  FixedList<cubicEqn, 4> hitEqn;
841  forAll(hitEqn, i)
842  {
843  hitEqn[i] = cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
844  }
845 
846  if (debug)
847  {
848  for (label i = 0; i < 4; ++ i)
849  {
850  Info<< (i ? " " : "Hit equation ") << i << " = "
851  << hitEqn[i] << endl;
852  }
853  Info<< " DetA equation = " << detA << endl;
854  }
855 
856  // Calculate the hit fraction
857  label iH = -1;
858  scalar muH = detA[0] > vSmall ? 1/detA[0] : vGreat;
859  for (label i = 0; i < 4; ++ i)
860  {
861  const Roots<3> mu = hitEqn[i].roots();
862 
863  for (label j = 0; j < 3; ++ j)
864  {
865  if
866  (
867  mu.type(j) == rootType::real
868  && hitEqn[i].derivative(mu[j]) < - vSmall
869  && hitEqn[i].derivative(mu[j]) < - mag(detA[0])*small
870  )
871  {
872  if (debug)
873  {
874  const barycentric yH
875  (
876  hitEqn[0].value(mu[j]),
877  hitEqn[1].value(mu[j]),
878  hitEqn[2].value(mu[j]),
879  hitEqn[3].value(mu[j])
880  );
881  const scalar detAH = detAEqn.value(mu[j]);
882 
883  Info<< "Hit on tet face " << i << " at local coordinate "
884  << (mag(detAH) > vSmall ? name(yH/detAH) : "???")
885  << ", " << mu[j]*detA[0]*100 << "% of the "
886  << "way along the track" << endl;
887  }
888 
889  if (0 <= mu[j] && mu[j] < muH)
890  {
891  iH = i;
892  muH = mu[j];
893  }
894  }
895  }
896  }
897 
898  // If there has been no hit on a degenerate or inverted tet then the
899  // displacement must be within the round off error. Advance the step
900  // fraction without moving and return.
901  if (iH == -1 && muH == vGreat)
902  {
903  stepFraction_ += fraction;
904  return 0;
905  }
906 
907  // Set the new coordinates
908  barycentric yH
909  (
910  hitEqn[0].value(muH),
911  hitEqn[1].value(muH),
912  hitEqn[2].value(muH),
913  hitEqn[3].value(muH)
914  );
915  // !!! <-- This fails if the tet collapses onto the particle, as detA tends
916  // to zero at the hit. In this instance, we can differentiate the hit and
917  // detA polynomials to find a limiting location, but this will not be on a
918  // triangle. We will then need to do a second track through the degenerate
919  // tet to find the final hit position. This second track is over zero
920  // distance and therefore can be of the static mesh type. This has not yet
921  // been implemented.
922  const scalar detAH = detAEqn.value(muH);
923  if (mag(detAH) < vSmall)
924  {
926  << "A moving tet collapsed onto a particle. This is not supported. "
927  << "The mesh is too poor, or the motion too severe, for particle "
928  << "tracking to function." << exit(FatalError);
929  }
930  yH /= detAH;
931 
932  // Clamp to zero any negative coordinates generated by round-off error
933  for (label i = 0; i < 4; ++ i)
934  {
935  yH.replace(i, i == iH ? 0 : max(0, yH[i]));
936  }
937 
938  // Re-normalise if within the tet
939  if (iH == -1)
940  {
941  yH /= cmptSum(yH);
942  }
943 
944  // Set the new position and hit index
945  coordinates_ = yH;
946  tetTriI = iH;
947 
948  // Set the proportion of the track that has been completed
949  stepFraction_ += fraction*muH*detA[0];
950 
951  if (debug)
952  {
953  if (iH != -1)
954  {
955  Info<< "Track hit tet face " << iH << " first" << endl;
956  }
957  else
958  {
959  Info<< "Track hit no tet faces" << endl;
960  }
961  Info<< "End local coordinates = " << yH << endl
962  << "End global coordinates = " << position() << endl
963  << "Tracking displacement = " << position() - x0 << endl
964  << muH*detA[0]*100 << "% of the step from "
965  << stepFraction_ - fraction*muH*detA[0] << " to "
966  << stepFraction_ - fraction*muH*detA[0] + fraction
967  << " completed" << endl << endl;
968  }
969 
970  // Accumulate fraction behind
971  if (muH*detA[0] < small || nTracksBehind_ > 0)
972  {
973  stepFractionBehind_ += (fraction != 0 ? fraction : 1)*muH*detA[0];
974 
975  if (stepFractionBehind_ > rootSmall)
976  {
977  stepFractionBehind_ = 0;
978  nTracksBehind_ = 0;
979  }
980  else
981  {
982  ++ nTracksBehind_;
983  }
984  }
985 
986  return iH != -1 ? 1 - muH*detA[0] : 0;
987 }
988 
989 
990 Foam::scalar Foam::particle::trackToTri
991 (
992  const vector& displacement,
993  const scalar fraction,
994  label& tetTriI
995 )
996 {
997  if (mesh_.moving() && (stepFraction_ != 1 || fraction != 0))
998  {
999  return trackToMovingTri(displacement, fraction, tetTriI);
1000  }
1001  else
1002  {
1003  return trackToStationaryTri(displacement, fraction, tetTriI);
1004  }
1005 }
1006 
1007 
1009 {
1010  if (cmptMin(mesh_.geometricD()) == -1)
1011  {
1012  vector pos = position(), posC = pos;
1013  meshTools::constrainToMeshCentre(mesh_, posC);
1014  return pos - posC;
1015  }
1016  else
1017  {
1018  return vector::zero;
1019  }
1020 }
1021 
1022 
1024 {}
1025 
1026 
1029  trackingData& td
1030 )
1031 {
1032  if (td.sendFromPatch == patch())
1033  {
1035  }
1036  else
1037  {
1039  (
1040  td.sendFromPatch,
1041  td.sendToPatchFace
1042  );
1043  }
1044 }
1045 
1046 
1048 {
1049  const polyPatch& pp = mesh_.boundaryMesh()[td.sendToPatch];
1050 
1051  if (isA<processorPolyPatch>(pp))
1052  {
1054  }
1055  else if (isA<nonConformalCyclicPolyPatch>(pp))
1056  {
1058  }
1059  else
1060  {
1062  << "Transfer patch type not recognised"
1063  << exit(FatalError);
1064  }
1065 }
1066 
1067 
1069 {
1070  // Store the local patch face in the face index
1071  facei_ = td.sendToPatchFace;
1072 }
1073 
1074 
1076 {
1077  const processorPolyPatch& ppp =
1078  refCast<const processorPolyPatch>(mesh_.boundaryMesh()[td.sendToPatch]);
1079 
1080  if (ppp.transform().transformsPosition())
1081  {
1083  }
1084 
1085  // Set the topology
1086  celli_ = ppp.faceCells()[facei_];
1087  facei_ += ppp.start();
1088  tetFacei_ = facei_;
1089 
1090  // Faces either side of a coupled patch are numbered in opposite
1091  // directions as their normals both point away from their connected
1092  // cells. The tet point therefore counts in the opposite direction from
1093  // the base point.
1094  tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
1095 
1096  // Reflect to account for the change of tri orientation in the new cell
1097  reflect();
1098 
1099  // Note that the position does not need transforming explicitly. The
1100  // face-triangle on the receive patch is the transformation of the one
1101  // on the send patch, so whilst the barycentric coordinates remain the
1102  // same, the change of triangle implicitly transforms the position.
1103 }
1104 
1105 
1108  const label sendFromPatch,
1109  const label sendToPatchFace
1110 )
1111 {
1112  const nonConformalCyclicPolyPatch& nccpp =
1113  static_cast<const nonConformalCyclicPolyPatch&>
1114  (mesh_.boundaryMesh()[sendFromPatch]);
1115 
1116  // Get the transformed position
1117  const vector pos = nccpp.transform().invTransformPosition(position());
1118 
1119  // Store the position in the barycentric data
1120  coordinates_ = barycentric(1 - cmptSum(pos), pos.x(), pos.y(), pos.z());
1121 
1122  // Break the topology
1123  celli_ = -1;
1124  tetFacei_ = -1;
1125  tetPti_ = -1;
1126 
1127  // Store the local patch face in the face index
1128  facei_ = sendToPatchFace;
1129 
1130  // Transform the properties
1131  if (nccpp.transform().transformsPosition())
1132  {
1134  }
1135 }
1136 
1137 
1140  const label sendToPatch
1141 )
1142 {
1143  const nonConformalCyclicPolyPatch& nccpp =
1144  static_cast<const nonConformalCyclicPolyPatch&>
1145  (mesh_.boundaryMesh()[sendToPatch]);
1146 
1147  // Get the position from the barycentric data
1148  const vector receivePos
1149  (
1150  coordinates_.b(),
1151  coordinates_.c(),
1152  coordinates_.d()
1153  );
1154 
1155  // Locate the particle on the receiving side
1156  locate
1157  (
1158  receivePos,
1159  mesh_.faceOwner()[facei_ + nccpp.origPatch().start()],
1160  false,
1161  "Particle crossed between " + nonConformalCyclicPolyPatch::typeName +
1162  " patches " + nccpp.name() + " and " + nccpp.nbrPatch().name() +
1163  " to a location outside of the mesh."
1164  );
1165 
1166  // The particle must remain associated with a face for the tracking to
1167  // register as incomplete
1168  facei_ = tetFacei_;
1169 }
1170 
1171 
1174  const transformer& transform
1175 )
1176 {
1177  // Get the transformed position
1178  const vector pos = transform.invTransformPosition(position());
1179 
1180  // Break the topology
1181  celli_ = -1;
1182  tetFacei_ = -1;
1183  tetPti_ = -1;
1184 
1185  // Store the position in the barycentric data
1186  coordinates_ = barycentric(1 - cmptSum(pos), pos.x(), pos.y(), pos.z());
1187 
1188  // Transform the properties
1189  if (transform.transformsPosition())
1190  {
1191  transformProperties(inv(transform));
1192  }
1193 }
1194 
1195 
1197 {
1198  // Get the position from the barycentric data
1199  const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
1200 
1201  // Create some arbitrary topology for the supplied cell
1202  celli_ = celli;
1203  tetFacei_ = mesh_.cells()[celli_][0];
1204  tetPti_ = 1;
1205 
1206  // Get the reverse transform and directly set the coordinates from the
1207  // position. This isn't likely to be correct; the particle is probably not
1208  // in this tet. It will, however, generate the correct vector when the
1209  // position method is called. A referred particle should never be tracked,
1210  // so this approximate topology is good enough. By using the nearby cell we
1211  // minimise the error associated with the incorrect topology.
1212  coordinates_ = barycentric(1, 0, 0, 0);
1213  if (mesh_.moving() && stepFraction_ != 1)
1214  {
1215  Pair<vector> centre;
1216  FixedList<scalar, 4> detA;
1218  movingTetReverseTransform(0, centre, detA, T);
1219  coordinates_ += (pos - centre[0]) & T[0]/detA[0];
1220  }
1221  else
1222  {
1223  vector centre;
1224  scalar detA;
1226  stationaryTetReverseTransform(centre, detA, T);
1227  coordinates_ += (pos - centre) & T/detA;
1228  }
1229 }
1230 
1231 
1234  const polyMesh& procMesh,
1235  const label procCell,
1236  const label procTetFace
1237 ) const
1238 {
1239  // The tet point on the procMesh differs from the current tet point if the
1240  // mesh and procMesh faces are of differing orientation. The change is the
1241  // same as in particle::correctAfterParallelTransfer.
1242 
1243  if
1244  (
1245  (mesh_.faceOwner()[tetFacei_] == celli_)
1246  == (procMesh.faceOwner()[procTetFace] == procCell)
1247  )
1248  {
1249  return tetPti_;
1250  }
1251  else
1252  {
1253  return procMesh.faces()[procTetFace].size() - 1 - tetPti_;
1254  }
1255 }
1256 
1257 
1260  const vector& position,
1261  const polyTopoChangeMap& mapper
1262 )
1263 {
1264  locate
1265  (
1266  position,
1267  mapper.reverseCellMap()[celli_],
1268  true,
1269  "Particle mapped to a location outside of the mesh."
1270  );
1271 }
1272 
1273 
1274 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
1275 
1276 bool Foam::operator==(const particle& pA, const particle& pB)
1277 {
1278  return (pA.origProc() == pB.origProc() && pA.origId() == pB.origId());
1279 }
1280 
1281 
1282 bool Foam::operator!=(const particle& pA, const particle& pB)
1283 {
1284  return !(pA == pB);
1285 }
1286 
1287 
1288 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
label origProc() const
Return the originating processor ID.
Definition: particleI.H:173
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
const word & name() const
Return name.
Vector-tensor class used to perform translations, rotations and scaling operations in 3D space...
Definition: transformer.H:83
vector deviationFromMeshCentre() const
Get the displacement from the mesh centre. Used to correct the.
Definition: particle.C:1008
bool moving() const
Is mesh moving.
Definition: polyMesh.H:525
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
void prepareForProcessorTransfer(trackingData &td)
Make changes prior to a transfer across a processor boundary.
Definition: particle.C:1068
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:54
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
label procTetPt(const polyMesh &procMesh, const label procCell, const label procTetFace) const
Return the tet point appropriate for decomposition or reconstruction.
Definition: particle.C:1233
void correctAfterParallelTransfer(trackingData &td)
Make changes following a parallel transfer. Runs either processor.
Definition: particle.C:1047
void correctAfterNonConformalCyclicTransfer(const label sendToPatch)
Make changes following a transfer across a non conformal cyclic.
Definition: particle.C:1139
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const Cmpt & b() const
Definition: BarycentricI.H:66
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
scalar trackToFace(const vector &displacement, const scalar fraction)
As particle::track, but stops when a face is hit.
Definition: particle.C:606
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int compare(const edge &, const edge &)
Compare edges.
Definition: edgeI.H:36
scalar trackToTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToFace, but stops when a tet triangle is hit. On.
Definition: particle.C:991
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
Cmpt cmptSum(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:500
void autoMap(const vector &position, const polyTopoChangeMap &mapper)
Map after a topology change.
Definition: particle.C:1259
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
dimensionedScalar y0(const dimensionedScalar &ds)
const Cmpt & z() const
Definition: VectorI.H:87
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
Definition: Roots.H:64
const cellList & cells() const
virtual const transformer & transform() const
Return transformation between the coupled patches.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
Definition: barycentric.H:45
void replace(const direction, const Cmpt &)
Definition: VectorSpaceI.H:149
scalar trackToCell(const vector &displacement, const scalar fraction)
As particle::track, but stops when a new cell is reached.
Definition: particle.C:584
label sendFromPatch
Patch from which to send the particle.
Definition: particle.H:114
virtual void transformProperties(const transformer &)
Transform the physical properties of the particle.
Definition: particle.C:1023
const dimensionedScalar c
Speed of light in a vacuum.
const Cmpt & y() const
Definition: VectorI.H:81
Base particle class.
Definition: particle.H:81
Neighbour processor patch.
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
Definition: polyMesh.C:894
dimensionedScalar pos(const dimensionedScalar &ds)
Cubic equation of the form a*x^3 + b*x^2 + c*x + d = 0.
Definition: cubicEqn.H:49
vector invTransformPosition(const vector &v) const
Inverse transform the given position.
Definition: transformerI.H:177
void prepareForInteractionListReferral(const transformer &transform)
Break the topology and store the cartesian position so that the.
Definition: particle.C:1173
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
void Swap(T &a, T &b)
Definition: Swap.H:43
const labelUList & faceCells() const
Return face-cell addressing.
Definition: polyPatch.C:340
tmp< fvMatrix< Type > > operator==(const fvMatrix< Type > &, const fvMatrix< Type > &)
void prepareForParallelTransfer(trackingData &td)
Make changes prior to a parallel transfer. Runs either processor or.
Definition: particle.C:1028
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
particle(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from components.
Definition: particle.C:474
3D tensor transformation operations.
face triFace(3)
void correctAfterInteractionListReferral(const label celli)
Correct the topology after referral. Locates the particle relative.
Definition: particle.C:1196
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1237
const Cmpt & d() const
Definition: BarycentricI.H:80
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1224
Non-conformal cyclic poly patch. As nonConformalCoupledPolyPatch, but the neighbouring patch is local...
const Cmpt & x() const
Definition: VectorI.H:75
label origId() const
Return the particle ID on the originating processor.
Definition: particleI.H:185
const Cmpt & c() const
Definition: BarycentricI.H:73
const polyPatch & origPatch() const
Original patch.
static const char nl
Definition: Ostream.H:260
const dimensionedScalar mu
Atomic mass unit.
bool onFace() const
Is the particle on a face?
Definition: particleI.H:254
#define FUNCTION_NAME
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
const labelList & reverseCellMap() const
Reverse cell map.
BarycentricTensor< scalar > barycentricTensor
A scalar version of the templated BarycentricTensor.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
bool transformsPosition() const
Return true if the transformer transforms a point.
Definition: transformerI.H:146
scalar trackToMovingTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToTri, but for moving meshes.
Definition: particle.C:794
virtual const transformer & transform() const
Return null transform between processor patches.
label sendToPatch
Patch to which to send the particle.
Definition: particle.H:117
void correctAfterProcessorTransfer(trackingData &td)
Make changes following a transfer across a processor boundary.
Definition: particle.C:1075
#define WarningInFunction
Report a warning using Foam::Warning.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:306
Templated 4x3 tensor derived from VectorSpace. Has 12 components. Can represent a barycentric transfo...
messageStream Info
const nonConformalCyclicPolyPatch & nbrPatch() const
Neighbour patch.
dimensioned< scalar > mag(const dimensioned< Type > &)
void prepareForNonConformalCyclicTransfer(const label sendToPatch, const label sendToPatchFace)
Make changes prior to a transfer across a non conformal cyclic.
Definition: particle.C:1107
bool onInternalFace() const
Is the particle on an internal face?
Definition: particleI.H:260
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
scalar trackToStationaryTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToTri, but for stationary meshes.
Definition: particle.C:659
volScalarField & p
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
scalar track(const vector &displacement, const scalar fraction)
Track along the displacement for a given fraction of the overall.
Definition: particle.C:560
bool operator!=(const particle &, const particle &)
Definition: particle.C:1282
scalar value(const scalar x) const
Evaluate at x.
Definition: cubicEqnI.H:103
void constrainToMeshCentre(const polyMesh &mesh, point &pt)
Set the constrained components of position to mesh centre.
Definition: meshTools.C:613
vector position() const
Return current particle position.
Definition: particleI.H:278
label patch() const
Return the index of patch that the particle is on.
Definition: particleI.H:272
Namespace for OpenFOAM.
void type(const direction i, const rootType t)
Set the type of the i-th root.
Definition: RootsI.H:120
label sendToPatchFace
Patch face to which to send the particle.
Definition: particle.H:120