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 std::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.edge(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);
430 <<
"Cell not found for particle position " << position <<
"."
437 const vector displacement =
439 - mesh_.cellCentres()[celli_]
440 +
vector(VSMALL, VSMALL, VSMALL);
445 const class cell&
c = mesh_.cells()[celli_];
446 scalar minF = VGREAT;
447 label minTetFacei = -1, minTetPti = -1;
450 const class face&
f = mesh_.faces()[
c[cellTetFacei]];
451 for (label tetPti = 1; tetPti <
f.size() - 1; ++tetPti)
454 tetFacei_ =
c[cellTetFacei];
459 const scalar
f = trackToTri(displacement, 0, tetTriI);
469 minTetFacei = tetFacei_;
478 tetFacei_ = minTetFacei;
481 track(displacement, 0);
494 static label nWarnings = 0;
495 static const label maxNWarnings = 100;
496 if ((nWarnings < maxNWarnings) && boundaryFail)
501 if (nWarnings == maxNWarnings)
504 <<
"Suppressing any further warnings about particles being "
505 <<
"located outside of the mesh." <<
endl;
520 const label tetFacei,
534 origId_(getNewParticleID())
546 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
555 origId_(getNewParticleID())
563 "Particle initialised with a location outside of the mesh"
573 const label tetFacei,
579 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
588 origId_(getNewParticleID())
598 "Particle initialised with a location outside of the mesh"
607 coordinates_(
p.coordinates_),
609 tetFacei_(
p.tetFacei_),
612 stepFraction_(
p.stepFraction_),
614 nBehind_(
p.nBehind_),
615 origProc_(
p.origProc_),
623 coordinates_(
p.coordinates_),
625 tetFacei_(
p.tetFacei_),
628 stepFraction_(
p.stepFraction_),
630 nBehind_(
p.nBehind_),
631 origProc_(
p.origProc_),
640 const vector& displacement,
641 const scalar fraction
644 scalar
f = trackToFace(displacement, fraction);
646 while (onInternalFace())
650 f *= trackToFace(
f*displacement,
f*fraction);
659 const vector& displacement,
660 const scalar fraction
665 label tetTriI = onFace() ? 0 : -1;
670 while (nBehind_ < maxNBehind_)
672 f *= trackToTri(
f*displacement,
f*fraction, tetTriI);
679 else if (tetTriI == 0)
693 static label stuckID = -1, stuckProc = -1;
694 if (origId_ != stuckID && origProc_ != stuckProc)
697 <<
"Particle #" << origId_ <<
" got stuck at " << position()
702 stuckProc = origProc_;
704 stepFraction_ +=
f*fraction;
715 const vector& displacement,
716 const scalar fraction,
720 const vector x0 = position();
721 const vector x1 = displacement;
726 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
727 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
734 stationaryTetReverseTransform(centre, detA,
T);
739 stationaryTetGeometry(o,
b, v1, v2);
740 Pout<<
"Tet points o=" << o <<
", b=" <<
b
741 <<
", v1=" << v1 <<
", v2=" << v2 <<
endl
742 <<
"Tet determinant = " << detA <<
endl
743 <<
"Start local coordinates = " <<
y0 <<
endl;
751 Pout<<
"Local displacement = " << Tx1 <<
"/" << detA <<
endl;
756 scalar muH = detA <= 0 ? VGREAT : 1/detA;
757 for (label i = 0; i < 4; ++ i)
759 if (Tx1[i] < - detA*SMALL)
761 scalar
mu = -
y0[i]/Tx1[i];
765 Pout<<
"Hit on tet face " << i <<
" at local coordinate "
766 <<
y0 +
mu*Tx1 <<
", " <<
mu*detA*100 <<
"% of the "
767 <<
"way along the track" <<
endl;
770 if (0 <=
mu &&
mu < muH)
782 for (label i = 0; i < 4; ++ i)
801 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
805 Pout<<
"Track hit no tet faces" <<
endl;
807 Pout<<
"End local coordinates = " << yH <<
endl
808 <<
"End global coordinates = " << position() <<
endl
809 <<
"Tracking displacement = " << position() - x0 <<
endl
810 << muH*detA*100 <<
"% of the step from " << stepFraction_ <<
" to "
811 << stepFraction_ + fraction <<
" completed" <<
endl <<
endl;
815 stepFraction_ += fraction*muH*detA;
819 Pout <<
"Step Fraction : " << stepFraction_ <<
endl;
820 Pout <<
"fraction*muH*detA : " << fraction*muH*detA <<
endl;
821 Pout <<
"static muH : " << muH <<
endl;
822 Pout <<
"origId() : " << origId() <<
endl;
826 if (detA <= 0 || nBehind_ > 0)
828 behind_ += muH*detA*
mag(displacement);
841 return iH != -1 ? 1 - muH*detA : 0;
847 const vector& displacement,
848 const scalar fraction,
852 const vector x0 = position();
853 const vector x1 = displacement;
858 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
859 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
866 movingTetReverseTransform(fraction, centre, detA,
T);
871 movingTetGeometry(fraction, o,
b, v1, v2);
872 Pout<<
"Tet points o=" << o[0] <<
", b=" <<
b[0]
873 <<
", v1=" << v1[0] <<
", v2=" << v2[0] <<
endl
874 <<
"Tet determinant = " << detA[0] <<
endl
875 <<
"Start local coordinates = " <<
y0[0] <<
endl;
880 const vector x0Rel = x0 - centre[0];
881 const vector x1Rel = x1 - centre[1];
884 cubicEqn detAEqn(
sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
887 ((x1Rel &
T[2]) + detA[3]*yC)*
sqr(detA[0]);
889 ((x1Rel &
T[1]) + (x0Rel &
T[2]) + detA[2]*yC)*detA[0];
891 ((x1Rel &
T[0]) + (x0Rel &
T[1]) + detA[1]*yC);
896 hitEqn[i] =
cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
901 for (label i = 0; i < 4; ++ i)
903 Pout<< (i ?
" " :
"Hit equation ") << i <<
" = "
904 << hitEqn[i] <<
endl;
906 Pout<<
" DetA equation = " << detA <<
endl;
911 scalar muH = detA[0] <= 0 ? VGREAT : 1/detA[0];
912 for (label i = 0; i < 4; ++i)
916 for (label j = 0; j < 3; ++j)
921 && hitEqn[i].derivative(
mu[j]) < - detA[0]*SMALL
928 hitEqn[0].value(
mu[j]),
929 hitEqn[1].value(
mu[j]),
930 hitEqn[2].value(
mu[j]),
931 hitEqn[3].value(
mu[j])
933 const scalar detAH = detAEqn.
value(
mu[j]);
935 Pout<<
"Hit on tet face " << i <<
" at local coordinate "
936 << (std::isnormal(detAH) ?
name(yH/detAH) :
"???")
937 <<
", " <<
mu[j]*detA[0]*100 <<
"% of the "
938 <<
"way along the track" <<
endl;
940 Pout<<
"derivative : " << hitEqn[i].derivative(
mu[j]) <<
nl
941 <<
" coord " << j <<
" mu[j]: " <<
mu[j] <<
nl
942 <<
" hitEq " << i <<
endl;
945 if (0 <=
mu[j] &&
mu[j] < muH)
957 hitEqn[0].value(muH),
958 hitEqn[1].value(muH),
959 hitEqn[2].value(muH),
969 const scalar detAH = detAEqn.
value(muH);
970 if (!std::isnormal(detAH))
973 <<
"A moving tet collapsed onto a particle. This is not supported. "
974 <<
"The mesh is too poor, or the motion too severe, for particle "
980 for (label i = 0; i < 4; ++ i)
995 scalar advance = muH*detA[0];
998 stepFraction_ += fraction*advance;
1001 if (detA[0] <= 0 || nBehind_ > 0)
1003 behind_ += muH*detA[0]*
mag(displacement);
1020 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
1024 Pout<<
"Track hit no tet faces" <<
endl;
1035 return iH != -1 ? 1 - muH*detA[0] : 0;
1041 const vector& displacement,
1042 const scalar fraction,
1046 if ((mesh_.moving() && (stepFraction_ != 1 || fraction != 0)))
1048 return trackToMovingTri(displacement, fraction, tetTriI);
1052 return trackToStationaryTri(displacement, fraction, tetTriI);
1059 if (
cmptMin(mesh_.geometricD()) == -1)
1083 facei_ = mesh_.boundaryMesh()[
patch()].whichFace(facei_);
1094 refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]);
1104 transformProperties(
T);
1114 transformProperties(-
s);
1119 facei_ += ppp.
start();
1124 tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
1165 const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
1169 tetFacei_ = mesh_.cells()[celli_][0];
1180 if (mesh_.moving() && stepFraction_ != 1)
1185 movingTetReverseTransform(0, centre, detA,
T);
1186 coordinates_ += (
pos - centre[0]) &
T[0]/detA[0];
1193 stationaryTetReverseTransform(centre, detA,
T);
1194 coordinates_ += (
pos - centre) &
T/detA;
1202 const label procCell,
1203 const label procTetFace
1212 (mesh_.faceOwner()[tetFacei_] == celli_)
1213 == (procMesh.
faceOwner()[procTetFace] == procCell)
1220 return procMesh.
faces()[procTetFace].size() - 1 - tetPti_;
1237 "Particle mapped to a location outside of the mesh"
1250 "Particle mapped to a location outside of the mesh"