68 const face&
f = mesh_.faces()[faceI];
71 label firstFullySubmergedPoint = -1;
76 scalar value = (mesh_.points()[
f[i]] - base) & normal;
77 if (
mag(value) < 10 * SMALL)
81 pointStatus_.append(value);
82 if (pointStatus_[i] > 10 * SMALL)
85 if (firstFullySubmergedPoint == -1)
87 firstFullySubmergedPoint = i;
92 if (inLiquid ==
f.size())
95 subFaceCentre_ = mesh_.faceCentres()[faceI];
96 subFaceArea_ = mesh_.faceAreas()[faceI];
99 else if (inLiquid == 0)
102 subFaceCentre_ =
Zero;
111 firstFullySubmergedPoint,
128 const scalar cutValue
134 label firstFullySubmergedPoint = -1;
140 pointStatus[i] = val[
f[i]] - cutValue;
141 if (
mag(pointStatus[i]) < 10 * SMALL)
145 if (pointStatus[i] > 10 * SMALL)
148 if (firstFullySubmergedPoint == -1)
150 firstFullySubmergedPoint = i;
155 if (inLiquid ==
f.size())
158 subFaceCentre_ =
f.centre(
points);
159 subFaceArea_ =
f.areaNormal(
points);
162 else if (inLiquid == 0)
165 subFaceCentre_ =
Zero;
175 firstFullySubmergedPoint,
227 const face&
f = mesh_.faces()[faceI];
230 if (
mag(Un0) > 1
e-12)
235 for (
const scalar fi :
f)
237 scalar value = ((mesh_.points()[fi] - x0) & n0) / Un0;
238 if (
mag(value) < 10 * SMALL)
242 pTimes_.append(value);
251 const label oldEdgeSign =
253 const label newEdgeSign =
256 if (newEdgeSign != oldEdgeSign)
262 if (nShifts == 2 || nShifts == 0)
264 dVf =
phi / magSf * timeIntegratedArea(faceI, dt, magSf, Un0);
266 else if (nShifts > 2)
271 fPts_tri[0] = mesh_.faceCentres()[faceI];
272 pTimes_tri[0] = ((fPts_tri[0] - x0) & n0)/Un0;
276 fPts_tri[1] = fPts[
pi];
277 pTimes_tri[1] = pTimes_[
pi];
279 pTimes_tri[2] = pTimes_[(
pi + 1) %
nPoints];
280 const scalar magSf_tri =
284 *(fPts_tri[2] - fPts_tri[0])
285 ^(fPts_tri[1] - fPts_tri[0])
288 const scalar phi_tri =
phi*magSf_tri/magSf;
306 <<
"Warning: nShifts = " << nShifts <<
" on face "
307 << faceI <<
" with pTimes = " << pTimes_
308 <<
" owned by cell " << mesh_.faceOwner()[faceI]
318 calcSubFace(faceI, -n0, x0);
319 const scalar alphaf =
mag(subFaceArea() / magSf);
324 <<
"Un0 is almost zero (" << Un0
325 <<
") - calculating dVf on face " << faceI
326 <<
" using subFaceFraction giving alphaf = " << alphaf
330 return phi * dt * alphaf;
353 pTimes_.append(times);
361 const label oldEdgeSign =
363 const label newEdgeSign =
366 if (newEdgeSign != oldEdgeSign)
374 dVf =
phi/magSf*timeIntegratedArea(faceI, dt, magSf, Un0);
392 scalar tIntArea = 0.0;
396 const scalar firstTime = pTimes[order.first()];
397 const scalar lastTime = pTimes[order.last()];
405 tIntArea = magSf * dt *
pos0(Un0);
416 tIntArea = magSf * dt * (1 -
pos0(Un0));
429 scalar initialArea = 0.0;
441 initialArea = magSf * (1.0 -
pos0(Un0));
442 tIntArea = initialArea * time;
443 cutPoints(fPts, pTimes, time, FIIL);
453 calcSubFace(
face(
identity(pTimes.size())), fPts, pTimes, time);
454 initialArea =
mag(subFaceArea());
455 cutPoints(fPts, pTimes, time, FIIL);
462 scalar prevTime = time;
463 const scalar tSmall =
max(1
e-6*dt, 10*SMALL);
465 for (
const scalar timeI : order)
467 if (timeI > prevTime + tSmall && timeI <= dt)
469 sortedTimes.append(timeI);
476 for (
const scalar newTime : sortedTimes)
480 cutPoints(fPts, pTimes, newTime, newFIIL);
484 quadAreaCoeffs(FIIL, newFIIL,
alpha,
beta);
486 tIntArea += (newTime - time) *
501 cutPoints(fPts, pTimes, dt, newFIIL);
505 quadAreaCoeffs(FIIL, newFIIL,
alpha,
beta);
507 tIntArea += (dt - time) *
514 tIntArea += magSf*(dt - lastTime)*
pos0(Un0);
521 void Foam::cutFaceAdvect::isoFacesToFile
530 fileName outputFile(filDir/(filNam +
".vtk"));
533 Info<<
"Writing file: " << outputFile <<
endl;
536 os <<
"# vtk DataFile Version 2.0" <<
nl
539 <<
"DATASET POLYDATA" <<
nl;
542 for (
const List<point>&
f : faces)
547 os <<
"POINTS " <<
nPoints <<
" float" <<
nl;
548 for (
const List<point>&
f : faces)
552 os <<
p.x() <<
' ' <<
p.y() <<
' ' <<
p.z() <<
nl;
557 << faces.size() <<
' ' << (
nPoints + faces.size()) <<
nl;
560 for (
const List<point>&
f : faces)
562 label endp =
f.size();
567 while (pointi < endp)
581 const scalar cutValue
585 const face&
f = mesh_.faces()[faceI];
587 label firstFullySubmergedPoint = -1;
592 scalar value = (
sign * pTimes_[i] - cutValue);
594 if (
mag(value) < 10 * SMALL)
598 pointStatus_.append(value);
599 if (pointStatus_[i] > 10 * SMALL)
602 if (firstFullySubmergedPoint == -1)
604 firstFullySubmergedPoint = i;
609 if (inLiquid ==
f.size())
612 subFaceCentre_ = mesh_.faceCentres()[faceI];
613 subFaceArea_ = mesh_.faceAreas()[faceI];
616 else if (inLiquid == 0)
619 subFaceCentre_ =
Zero;
628 firstFullySubmergedPoint,
649 scalar tIntArea = 0.0;
653 const scalar firstTime = pTimes_[order.first()];
654 const scalar lastTime = pTimes_[order.last()];
662 tIntArea = magSf* dt *
pos0(Un0);
673 tIntArea = magSf * dt * (1 -
pos0(Un0));
686 scalar initialArea = 0.0;
698 initialArea = magSf * (1.0 -
pos0(Un0));
699 tIntArea = initialArea * time;
700 cutPoints(faceI, time, FIIL);
708 calcSubFace(faceI, -
sign(Un0), time);
709 initialArea =
mag(subFaceArea());
710 cutPoints(faceI, time, FIIL);
717 scalar prevTime = time;
718 const scalar tSmall =
max(1
e-6*dt, 10*SMALL);
719 for (
const label oI : order)
721 const scalar timeI = pTimes_[oI];
722 if (timeI > prevTime + tSmall && timeI <= dt)
724 sortedTimes.append(timeI);
731 for (
const scalar newTime : sortedTimes)
735 cutPoints(faceI, newTime, newFIIL);
740 quadAreaCoeffs(FIIL, newFIIL,
alpha,
beta);
744 * (initialArea +
sign(Un0)
759 cutPoints(faceI, dt, newFIIL);
763 quadAreaCoeffs(FIIL, newFIIL,
alpha,
beta);
773 tIntArea += magSf * (dt - lastTime) *
pos0(Un0);
789 const label np0 = pf0.size();
790 const label np1 = pf1.size();
826 if (((
B -
A) & (
D -
C)) > 0)
837 const scalar Bx =
mag(
B -
A);
845 else if (
mag(
C -
D) > 10 * SMALL)
857 yhat -= ((yhat & xhat) * xhat);
859 if (
mag(yhat) > 10 * SMALL)
863 const scalar Cx = (
C -
A) & xhat;
864 const scalar Cy =
mag((
C -
A) & yhat);
865 const scalar Dx = (
D -
A) & xhat;
866 const scalar Dy =
mag((
D -
A) & yhat);
869 alpha = 0.5 * ((Cx - Bx) * Dy - Dx * Cy);
870 beta = 0.5 * Bx * (Dy + Cy);
876 <<
"Vertex face was cut at " << pf0 <<
" and at "
889 const face&
f = mesh_.faces()[faceI];
891 scalar f1(pTimes_[0]);
894 if (
mag(f1 - f0) < 10 * SMALL)
902 scalar f2 = pTimes_[pi2];
905 if (
mag(f2 - f0) < 10 * SMALL)
910 if ((f1 < f0 && f2 > f0) || (f1 > f0 && f2 < f0))
912 const scalar
s = (f0 - f1) / (f2 - f1);
915 mesh_.points()[
f[
pi]]
916 +
s*(mesh_.points()[
f[pi2]] - mesh_.points()[
f[
pi]])
926 if (cutPoints.size() > 2)
929 <<
"cutPoints = " << cutPoints
930 <<
" for pts = " <<
f.points(mesh_.points())
931 <<
", f - f0 = " <<
f - f0 <<
" and f0 = " << f0
945 const label
nPoints = pts.size();
949 if (
mag(f1 - f0) < 10 * SMALL)
960 if (
mag(f2 - f0) < 10 * SMALL)
965 if ((f1 < f0 && f2 > f0) || (f1 > f0 && f2 < f0))
967 const scalar
s = (f0 - f1) / (f2 - f1);
968 cutPoints.
append(pts[
pi] +
s * (pts[pi2] - pts[
pi]));
977 if (cutPoints.size() > 2)
980 <<
"cutPoints = " << cutPoints <<
" for pts = " << pts
981 <<
", f - f0 = " <<
f - f0 <<
" and f0 = " << f0
989 return subFaceCentre_;
1001 return subFacePoints_;
1007 return surfacePoints_;
1013 subFaceCentre_ =
Zero;
1014 subFaceArea_ =
Zero;
1015 subFacePoints_.clear();
1016 surfacePoints_.clear();
1017 pointStatus_.clear();