43 const Foam::label Foam::particle::maxNBehind_ = 10;
56 "writeLagrangianPositions",
64 void Foam::particle::stationaryTetReverseTransform
81 detA = ab & (ac ^ ad);
93 void Foam::particle::movingTetReverseTransform
95 const scalar fraction,
97 FixedList<scalar, 4>& detA,
98 FixedList<barycentricTensor, 3>&
T
101 Pair<barycentricTensor>
A = movingTetTransform(fraction);
103 const Pair<vector> ab(
A[0].
b() -
A[0].a(),
A[1].
b() -
A[1].a());
104 const Pair<vector> ac(
A[0].
c() -
A[0].a(),
A[1].
c() -
A[1].a());
105 const Pair<vector> ad(
A[0].d() -
A[0].a(),
A[1].d() -
A[1].a());
106 const Pair<vector> bc(
A[0].
c() -
A[0].
b(),
A[1].
c() -
A[1].
b());
107 const Pair<vector> bd(
A[0].d() -
A[0].
b(),
A[1].d() -
A[1].
b());
109 centre[0] =
A[0].a();
110 centre[1] =
A[1].a();
112 detA[0] = ab[0] & (ac[0] ^ ad[0]);
114 (ab[1] & (ac[0] ^ ad[0]))
115 + (ab[0] & (ac[1] ^ ad[0]))
116 + (ab[0] & (ac[0] ^ ad[1]));
118 (ab[0] & (ac[1] ^ ad[1]))
119 + (ab[1] & (ac[0] ^ ad[1]))
120 + (ab[1] & (ac[1] ^ ad[0]));
121 detA[3] = ab[1] & (ac[1] ^ ad[1]);
132 (bd[0] ^ bc[1]) + (bd[1] ^ bc[0]),
133 (ac[0] ^ ad[1]) + (ac[1] ^ ad[0]),
134 (ad[0] ^ ab[1]) + (ad[1] ^ ab[0]),
135 (ab[0] ^ ac[1]) + (ab[1] ^ ac[0])
147 void Foam::particle::reflect()
149 Swap(coordinates_.c(), coordinates_.d());
153 void Foam::particle::rotate(
const bool reverse)
157 scalar temp = coordinates_.b();
158 coordinates_.b() = coordinates_.c();
159 coordinates_.c() = coordinates_.d();
160 coordinates_.d() = temp;
164 scalar temp = coordinates_.d();
165 coordinates_.d() = coordinates_.c();
166 coordinates_.c() = coordinates_.b();
167 coordinates_.b() = temp;
172 void Foam::particle::changeTet(
const label tetTriI)
174 const bool isOwner = mesh_.faceOwner()[tetFacei_] == celli_;
176 const label firstTetPtI = 1;
177 const label lastTetPtI = mesh_.faces()[tetFacei_].size() - 2;
183 else if (tetTriI == 2)
187 if (tetPti_ == lastTetPtI)
199 if (tetPti_ == firstTetPtI)
210 else if (tetTriI == 3)
214 if (tetPti_ == firstTetPtI)
226 if (tetPti_ == lastTetPtI)
240 <<
"Changing tet without changing cell should only happen when the "
241 <<
"track is on triangle 1, 2 or 3."
247 void Foam::particle::changeFace(
const label tetTriI)
250 const triFace triOldIs(currentTetIndices().faceTriIs(mesh_));
256 sharedEdge = edge(triOldIs[1], triOldIs[2]);
258 else if (tetTriI == 2)
260 sharedEdge = edge(triOldIs[2], triOldIs[0]);
262 else if (tetTriI == 3)
264 sharedEdge = edge(triOldIs[0], triOldIs[1]);
269 <<
"Changing face without changing cell should only happen when the"
270 <<
" track is on triangle 1, 2 or 3."
273 sharedEdge = edge(-1, -1);
279 forAll(mesh_.cells()[celli_], cellFaceI)
281 const label newFaceI = mesh_.cells()[celli_][cellFaceI];
282 const class face& newFace = mesh_.faces()[newFaceI];
283 const label newOwner = mesh_.faceOwner()[newFaceI];
286 if (tetFacei_ == newFaceI)
294 const label edgeComp = newOwner == celli_ ? -1 : +1;
299 edgeI < newFace.size()
300 &&
edge::compare(sharedEdge, newFace.faceEdge(edgeI)) != edgeComp;
305 if (edgeI >= newFace.size())
311 const label newBaseI =
max(0, mesh_.tetBasePtIs()[newFaceI]);
312 edgeI = (edgeI - newBaseI + newFace.size()) % newFace.size();
317 edgeI =
min(
max(1, edgeI), newFace.size() - 2);
320 tetFacei_ = newFaceI;
330 <<
"The search for an edge-connected face and tet-point failed."
335 if (sharedEdge.otherVertex(triOldIs[1]) == -1)
339 else if (sharedEdge.otherVertex(triOldIs[2]) == -1)
345 const triFace triNewIs = currentTetIndices().faceTriIs(mesh_);
351 if (sharedEdge.otherVertex(triNewIs[1]) == -1)
355 else if (sharedEdge.otherVertex(triNewIs[2]) == -1)
362 void Foam::particle::changeCell()
365 const label ownerCellI = mesh_.faceOwner()[tetFacei_];
366 const bool isOwner = celli_ == ownerCellI;
367 celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
373 void Foam::particle::changeToMasterPatch()
375 label thisPatch =
patch();
377 forAll(mesh_.cells()[celli_], cellFaceI)
380 const label otherFaceI = mesh_.cells()[celli_][cellFaceI];
381 if (facei_ == otherFaceI || mesh_.isInternalFace(otherFaceI))
389 const class face& thisFace = mesh_.faces()[facei_];
390 const class face&
otherFace = mesh_.faces()[otherFaceI];
393 const label otherPatch =
394 mesh_.boundaryMesh().whichPatch(otherFaceI);
395 if (thisPatch > otherPatch)
398 thisPatch = otherPatch;
407 void Foam::particle::locate
412 const bool boundaryFail,
413 const string& boundaryMsg
426 celli_ = mesh_.cellTree().findInside(position);
431 <<
"Cell not found for particle position " << position <<
"."
435 const vector displacement = position - mesh_.cellCentres()[celli_];
440 const class cell&
c = mesh_.cells()[celli_];
441 scalar minF = VGREAT;
442 label minTetFacei = -1, minTetPti = -1;
445 const class face&
f = mesh_.faces()[
c[cellTetFacei]];
446 for (label tetPti = 1; tetPti <
f.size() - 1; ++tetPti)
449 tetFacei_ =
c[cellTetFacei];
454 const scalar
f = trackToTri(displacement, 0, tetTriI);
464 minTetFacei = tetFacei_;
473 tetFacei_ = minTetFacei;
476 track(displacement, 0);
489 static label nWarnings = 0;
490 static const label maxNWarnings = 100;
491 if ((nWarnings < maxNWarnings) && boundaryFail)
496 if (nWarnings == maxNWarnings)
499 <<
"Suppressing any further warnings about particles being "
500 <<
"located outside of the mesh." <<
endl;
515 const label tetFacei,
529 origId_(getNewParticleID())
541 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
550 origId_(getNewParticleID())
558 "Particle initialised with a location outside of the mesh"
568 const label tetFacei,
574 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
583 origId_(getNewParticleID())
593 "Particle initialised with a location outside of the mesh"
602 coordinates_(
p.coordinates_),
604 tetFacei_(
p.tetFacei_),
607 stepFraction_(
p.stepFraction_),
609 nBehind_(
p.nBehind_),
610 origProc_(
p.origProc_),
618 coordinates_(
p.coordinates_),
620 tetFacei_(
p.tetFacei_),
623 stepFraction_(
p.stepFraction_),
625 nBehind_(
p.nBehind_),
626 origProc_(
p.origProc_),
635 const vector& displacement,
636 const scalar fraction
639 scalar
f = trackToFace(displacement, fraction);
641 while (onInternalFace())
645 f *= trackToFace(
f*displacement,
f*fraction);
654 const vector& displacement,
655 const scalar fraction
660 label tetTriI = onFace() ? 0 : -1;
665 while (nBehind_ < maxNBehind_)
667 f *= trackToTri(
f*displacement,
f*fraction, tetTriI);
674 else if (tetTriI == 0)
688 static label stuckID = -1, stuckProc = -1;
689 if (origId_ != stuckID && origProc_ != stuckProc)
692 <<
"Particle #" << origId_ <<
" got stuck at " << position()
697 stuckProc = origProc_;
699 stepFraction_ +=
f*fraction;
710 const vector& displacement,
711 const scalar fraction,
715 const vector x0 = position();
716 const vector x1 = displacement;
721 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
722 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
729 stationaryTetReverseTransform(centre, detA,
T);
734 stationaryTetGeometry(o,
b, v1, v2);
735 Pout<<
"Tet points o=" << o <<
", b=" <<
b
736 <<
", v1=" << v1 <<
", v2=" << v2 <<
endl
737 <<
"Tet determinant = " << detA <<
endl
738 <<
"Start local coordinates = " <<
y0 <<
endl;
746 Pout<<
"Local displacement = " << Tx1 <<
"/" << detA <<
endl;
751 scalar muH = std::isnormal(detA) && detA <= 0 ? VGREAT : 1/detA;
752 for (label i = 0; i < 4; ++ i)
754 if (Tx1[i] < - detA*SMALL)
756 scalar
mu = -
y0[i]/Tx1[i];
760 Pout<<
"Hit on tet face " << i <<
" at local coordinate "
761 <<
y0 +
mu*Tx1 <<
", " <<
mu*detA*100 <<
"% of the "
762 <<
"way along the track" <<
endl;
765 if (0 <=
mu &&
mu < muH)
777 for (label i = 0; i < 4; ++ i)
796 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
800 Pout<<
"Track hit no tet faces" <<
endl;
802 Pout<<
"End local coordinates = " << yH <<
endl
803 <<
"End global coordinates = " << position() <<
endl
804 <<
"Tracking displacement = " << position() - x0 <<
endl
805 << muH*detA*100 <<
"% of the step from " << stepFraction_ <<
" to "
806 << stepFraction_ + fraction <<
" completed" <<
endl <<
endl;
810 stepFraction_ += fraction*muH*detA;
814 Info <<
"Step Fraction : " << stepFraction_ <<
endl;
815 Info <<
"fraction*muH*detA : " << fraction*muH*detA <<
endl;
816 Info <<
"static muH : " << muH <<
endl;
817 Info <<
"origId() : " << origId() <<
endl;
821 if (detA <= 0 || nBehind_ > 0)
823 behind_ += muH*detA*
mag(displacement);
836 return iH != -1 ? 1 - muH*detA : 0;
842 const vector& displacement,
843 const scalar fraction,
847 const vector x0 = position();
848 const vector x1 = displacement;
853 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
854 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
861 movingTetReverseTransform(fraction, centre, detA,
T);
866 movingTetGeometry(fraction, o,
b, v1, v2);
867 Info<<
"Tet points o=" << o[0] <<
", b=" <<
b[0]
868 <<
", v1=" << v1[0] <<
", v2=" << v2[0] <<
endl
869 <<
"Tet determinant = " << detA[0] <<
endl
870 <<
"Start local coordinates = " <<
y0[0] <<
endl;
875 const vector x0Rel = x0 - centre[0];
876 const vector x1Rel = x1 - centre[1];
879 cubicEqn detAEqn(
sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
882 ((x1Rel &
T[2]) + detA[3]*yC)*
sqr(detA[0]);
884 ((x1Rel &
T[1]) + (x0Rel &
T[2]) + detA[2]*yC)*detA[0];
886 ((x1Rel &
T[0]) + (x0Rel &
T[1]) + detA[1]*yC);
891 hitEqn[i] =
cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
896 for (label i = 0; i < 4; ++ i)
898 Info<< (i ?
" " :
"Hit equation ") << i <<
" = "
899 << hitEqn[i] <<
endl;
901 Info<<
" DetA equation = " << detA <<
endl;
906 scalar muH = std::isnormal(detA[0]) && detA[0] <= 0 ? VGREAT : 1/detA[0];
907 for (label i = 0; i < 4; ++i)
911 for (label j = 0; j < 3; ++j)
916 && hitEqn[i].derivative(
mu[j]) < - detA[0]*SMALL
923 hitEqn[0].value(
mu[j]),
924 hitEqn[1].value(
mu[j]),
925 hitEqn[2].value(
mu[j]),
926 hitEqn[3].value(
mu[j])
928 const scalar detAH = detAEqn.
value(
mu[j]);
930 Info<<
"Hit on tet face " << i <<
" at local coordinate "
931 << (std::isnormal(detAH) ?
name(yH/detAH) :
"???")
932 <<
", " <<
mu[j]*detA[0]*100 <<
"% of the "
933 <<
"way along the track" <<
endl;
935 Info<<
"derivative : " << hitEqn[i].derivative(
mu[j]) <<
nl
936 <<
" coord " << j <<
" mu[j]: " <<
mu[j] <<
nl
937 <<
" hitEq " << i <<
endl;
940 if (0 <=
mu[j] &&
mu[j] < muH)
952 hitEqn[0].value(muH),
953 hitEqn[1].value(muH),
954 hitEqn[2].value(muH),
964 const scalar detAH = detAEqn.
value(muH);
965 if (!std::isnormal(detAH))
968 <<
"A moving tet collapsed onto a particle. This is not supported. "
969 <<
"The mesh is too poor, or the motion too severe, for particle "
975 for (label i = 0; i < 4; ++ i)
990 scalar advance = muH*detA[0];
993 stepFraction_ += fraction*advance;
996 if (detA[0] <= 0 || nBehind_ > 0)
998 behind_ += muH*detA[0]*
mag(displacement);
1015 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
1019 Pout<<
"Track hit no tet faces" <<
endl;
1030 return iH != -1 ? 1 - muH*detA[0] : 0;
1036 const vector& displacement,
1037 const scalar fraction,
1041 if ((mesh_.moving() && (stepFraction_ != 1 || fraction != 0)))
1043 return trackToMovingTri(displacement, fraction, tetTriI);
1047 return trackToStationaryTri(displacement, fraction, tetTriI);
1054 if (
cmptMin(mesh_.geometricD()) == -1)
1078 facei_ = mesh_.boundaryMesh()[
patch()].whichFace(facei_);
1089 refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]);
1099 transformProperties(
T);
1109 transformProperties(-
s);
1114 facei_ += ppp.
start();
1119 tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
1160 const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
1164 tetFacei_ = mesh_.cells()[celli_][0];
1175 if (mesh_.moving() && stepFraction_ != 1)
1180 movingTetReverseTransform(0, centre, detA,
T);
1181 coordinates_ += (
pos - centre[0]) &
T[0]/detA[0];
1188 stationaryTetReverseTransform(centre, detA,
T);
1189 coordinates_ += (
pos - centre) &
T/detA;
1197 const label procCell,
1198 const label procTetFace
1207 (mesh_.faceOwner()[tetFacei_] == celli_)
1208 == (procMesh.
faceOwner()[procTetFace] == procCell)
1215 return procMesh.
faces()[procTetFace].size() - 1 - tetPti_;
1232 "Particle mapped to a location outside of the mesh"
1245 "Particle mapped to a location outside of the mesh"