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;
275 fPts_tri[1] = fPts[
pi];
276 pTimes_tri[1] = pTimes_[
pi];
278 pTimes_tri[2] = pTimes_[(
pi + 1) %
nPoints];
279 const scalar magSf_tri =
283 *(fPts_tri[2] - fPts_tri[0])
284 ^(fPts_tri[1] - fPts_tri[0])
287 const scalar phi_tri =
phi*magSf_tri/magSf;
305 <<
"Warning: nShifts = " << nShifts <<
" on face " 306 << faceI <<
" with pTimes = " << pTimes_
307 <<
" owned by cell " << mesh_.faceOwner()[faceI]
317 calcSubFace(faceI, -n0, x0);
318 const scalar alphaf =
mag(subFaceArea() / magSf);
323 <<
"Un0 is almost zero (" << Un0
324 <<
") - calculating dVf on face " << faceI
325 <<
" using subFaceFraction giving alphaf = " << alphaf
329 return phi * dt * alphaf;
352 pTimes_.append(times);
360 const label oldEdgeSign =
362 const label newEdgeSign =
365 if (newEdgeSign != oldEdgeSign)
373 dVf =
phi/magSf*timeIntegratedArea(faceI, dt, magSf, Un0);
391 scalar tIntArea = 0.0;
395 const scalar firstTime = pTimes[order.first()];
396 const scalar lastTime = pTimes[order.last()];
404 tIntArea = magSf * dt *
pos0(Un0);
415 tIntArea = magSf * dt *
neg(Un0);
426 DynamicList<point> FIIL(3);
428 scalar initialArea = 0.0;
440 initialArea = magSf *
neg(Un0);
441 tIntArea = initialArea * time;
442 cutPoints(fPts, pTimes, time, FIIL);
452 calcSubFace(face(
identity(pTimes.size())), fPts, pTimes, time);
453 initialArea =
mag(subFaceArea());
454 cutPoints(fPts, pTimes, time, FIIL);
459 DynamicList<scalar> sortedTimes(pTimes.size());
461 scalar prevTime = time;
462 const scalar tSmall =
max(1
e-6*dt, 10*SMALL);
464 for (
const scalar timeI : order)
466 if (timeI > prevTime + tSmall && timeI <= dt)
468 sortedTimes.append(timeI);
475 for (
const scalar newTime : sortedTimes)
478 DynamicList<point> newFIIL(3);
479 cutPoints(fPts, pTimes, newTime, newFIIL);
483 quadAreaCoeffs(FIIL, newFIIL,
alpha,
beta);
485 tIntArea += (newTime - time) *
499 DynamicList<point> newFIIL(3);
500 cutPoints(fPts, pTimes, dt, newFIIL);
504 quadAreaCoeffs(FIIL, newFIIL,
alpha,
beta);
506 tIntArea += (dt - time) *
513 tIntArea += magSf*(dt - lastTime)*
pos0(Un0);
520 void Foam::cutFaceAdvect::isoFacesToFile
522 const DynamicList<List<point>>& faces,
529 fileName outputFile(filDir/(filNam +
".vtk"));
532 Info<<
"Writing file: " << outputFile <<
endl;
534 OFstream
os(outputFile);
535 os <<
"# vtk DataFile Version 2.0" <<
nl 538 <<
"DATASET POLYDATA" <<
nl;
541 for (
const List<point>&
f : faces)
547 for (
const List<point>&
f : faces)
551 os <<
p.x() <<
' ' <<
p.y() <<
' ' <<
p.z() <<
nl;
556 << faces.size() <<
' ' << (
nPoints + faces.size()) <<
nl;
559 for (
const List<point>&
f : faces)
561 label endp =
f.
size();
566 while (pointi < endp)
580 const scalar cutValue
584 const face&
f = mesh_.faces()[faceI];
586 label firstFullySubmergedPoint = -1;
591 scalar value = (
sign * pTimes_[i] - cutValue);
593 if (
mag(value) < 10 * SMALL)
597 pointStatus_.append(value);
598 if (pointStatus_[i] > 10 * SMALL)
601 if (firstFullySubmergedPoint == -1)
603 firstFullySubmergedPoint = i;
608 if (inLiquid ==
f.
size())
611 subFaceCentre_ = mesh_.faceCentres()[faceI];
612 subFaceArea_ = mesh_.faceAreas()[faceI];
615 else if (inLiquid == 0)
618 subFaceCentre_ =
Zero;
627 firstFullySubmergedPoint,
648 scalar tIntArea = 0.0;
652 const scalar firstTime = pTimes_[order.first()];
653 const scalar lastTime = pTimes_[order.last()];
661 tIntArea = magSf* dt *
pos0(Un0);
672 tIntArea = magSf * dt *
neg(Un0);
683 DynamicList<point> FIIL(3);
685 scalar initialArea = 0.0;
697 initialArea = magSf *
neg(Un0);
698 tIntArea = initialArea * time;
699 cutPoints(faceI, time, FIIL);
707 calcSubFace(faceI, -
sign(Un0), time);
708 initialArea =
mag(subFaceArea());
709 cutPoints(faceI, time, FIIL);
714 DynamicList<scalar> sortedTimes(pTimes_.size());
716 scalar prevTime = time;
717 const scalar tSmall =
max(1
e-6*dt, 10*SMALL);
718 for (
const label oI : order)
720 const scalar timeI = pTimes_[oI];
721 if (timeI > prevTime + tSmall && timeI <= dt)
723 sortedTimes.append(timeI);
730 for (
const scalar newTime : sortedTimes)
733 DynamicList<point> newFIIL(3);
734 cutPoints(faceI, newTime, newFIIL);
739 quadAreaCoeffs(FIIL, newFIIL,
alpha,
beta);
743 * (initialArea +
sign(Un0)
757 DynamicList<point> newFIIL(3);
758 cutPoints(faceI, dt, newFIIL);
762 quadAreaCoeffs(FIIL, newFIIL,
alpha,
beta);
772 tIntArea += magSf * (dt - lastTime) *
pos0(Un0);
781 const DynamicList<point>& pf0,
782 const DynamicList<point>& pf1,
788 const label np0 = pf0.size();
789 const label np1 = pf1.size();
825 if (((
B -
A) & (
D -
C)) > 0)
836 const scalar Bx =
mag(
B -
A);
844 else if (
mag(
C -
D) > 10 * SMALL)
856 yhat -= ((yhat & xhat) * xhat);
858 if (
mag(yhat) > 10 * SMALL)
862 const scalar Cx = (
C -
A) & xhat;
863 const scalar Cy =
mag((
C -
A) & yhat);
864 const scalar Dx = (
D -
A) & xhat;
865 const scalar Dy =
mag((
D -
A) & yhat);
868 alpha = 0.5 * ((Cx - Bx) * Dy - Dx * Cy);
869 beta = 0.5 * Bx * (Dy + Cy);
875 <<
"Vertex face was cut at " << pf0 <<
" and at " 885 DynamicList<point>& cutPoints
888 const face&
f = mesh_.faces()[faceI];
890 scalar f1(pTimes_[0]);
893 if (
mag(f1 - f0) < 10 * SMALL)
901 scalar f2 = pTimes_[pi2];
904 if (
mag(f2 - f0) < 10 * SMALL)
909 if ((f1 < f0 && f2 > f0) || (f1 > f0 && f2 < f0))
911 const scalar
s = (f0 - f1) / (f2 - f1);
914 mesh_.points()[
f[
pi]]
915 +
s*(mesh_.points()[
f[pi2]] - mesh_.points()[
f[
pi]])
920 cutPoints.append(mesh_.points()[
f[
pi]]);
925 if (cutPoints.size() > 2)
928 <<
"cutPoints = " << cutPoints
929 <<
" for pts = " <<
f.points(mesh_.points())
930 <<
", f - f0 = " <<
f - f0 <<
" and f0 = " << f0
941 DynamicList<point>& cutPoints
948 if (
mag(f1 - f0) < 10 * SMALL)
959 if (
mag(f2 - f0) < 10 * SMALL)
964 if ((f1 < f0 && f2 > f0) || (f1 > f0 && f2 < f0))
966 const scalar
s = (f0 - f1) / (f2 - f1);
971 cutPoints.append(
pts[
pi]);
976 if (cutPoints.size() > 2)
979 <<
"cutPoints = " << cutPoints <<
" for pts = " <<
pts 980 <<
", f - f0 = " <<
f - f0 <<
" and f0 = " << f0
988 subFaceCentre_ =
Zero;
990 subFacePoints_.clear();
991 surfacePoints_.clear();
992 pointStatus_.clear();
dimensionedScalar sign(const dimensionedScalar &ds)
scalar timeIntegratedFaceFlux(const label faceI, const vector &x0, const vector &n0, const scalar Un0, const scalar dt, const scalar phi, const scalar magSf)
Calculate time integrated flux for a face.
void size(const label n)
Older name for setAddressableSize.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label calcSubFace(const label faceI, const vector &normal, const vector &base)
Calculates cut centre and cut area (plicReconstruction)
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar neg(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
Base class for cutting a face, faceI, of an fvMesh, mesh_, at its intersections.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
constexpr scalar pi(M_PI)
dimensionedScalar pos0(const dimensionedScalar &ds)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
OBJstream os(runTime.globalPath()/outputName)
void cutPoints(const label faceI, const scalar f0, DynamicList< point > &cutPoints)
Get cutPoints of face.
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
Mesh data needed to do the Finite Volume discretisation.
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
void quadAreaCoeffs(const DynamicPointList &pf0, const DynamicPointList &pf1, scalar &quadArea, scalar &intQuadArea) const
For face with vertices p calculate its area and integrated area.
const dimensionedScalar & D
List< label > labelList
A List of labels.
void calcSubFace(const label faceI, const scalarList &pointStatus, label firstFullySubmergedPoint, DynamicList< point > &subFacePoints, DynamicList< point > &surfacePoints, label &faceStatus, vector &subFaceCentre, vector &subFaceArea)
Calculate cut points along edges of face with pointStatus, pointfield and computes geometric informat...
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
cutFaceAdvect(const fvMesh &mesh, const volScalarField &alpha1)
Construct from fvMesh and a scalarField.
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
static constexpr const zero Zero
Global zero (0)
scalar timeIntegratedArea(const label faceI, const scalar dt, const scalar magSf, const scalar Un0)
Calculate time integrated area for a face.
const volScalarField & alpha1
void clearStorage()
Resets internal variables.