59 <<
"normal length is zero in cell: " << celli <<
nl 60 <<
"try increasing nCorrectors" <<
endl;
72 fvert[
pi] = (mesh_.
points()[pLabels[
pi]] - mesh_.
C()[celli]) & (normal);
76 scalar f1 = fvert[order.first()];
77 scalar f2 = fvert[order.last()];
92 label L2 = fvert.size() - 1;
99 L3 = round(0.5*(L1 + L2));
100 f3 = fvert[order[L3]];
105 L1 = L3; f1 = f3; a1 = a3;
109 L2 = L3; f2 = f3; a2 = a3;
117 if (
mag(f1 - f2) < 10*SMALL)
122 if (
mag(a1 - a2) < tol)
124 return cutCell_.
calcSubCell(celli, 0.5*(f1 + f2), normal);
133 f3 = f1 + (f2 - f1)/3;
137 scalar f4 = f1 + 2*(f2 - f1)/3;
144 a[0] = a1, a[1] = a3, a[2] = a4, a[3] = a2;
145 fOld[0] = f1, fOld[1] = f3, fOld[2] = f4, fOld[3] = f2;
146 f[0] = 0,
f[1] = (f3-f1)/(f2-f1),
f[2] = (f4-f1)/(f2-f1),
f[3] = 1;
152 M[i][j] =
pow(
f[i], 3 - j);
162 f3 =
f[1]; a3 = a[1];
165 while (res > tol && nIter < 10*maxIter)
168 /(3*
C[0]*
sqr(f3) + 2*
C[1]*f3 +
C[2]);
169 a3 =
C[0]*
pow3(f3) +
C[1]*
sqr(f3) +
C[2]*f3 +
C[3];
174 f3 = f3*(f2 - f1) + f1;
177 label status = cutCell_.
calcSubCell(celli, f3, normal);
194 scalar x1 =
max(1
e-3*(f2 - f1), 100*SMALL);
202 while (res > tol && nIter < maxIter && g1 != g2)
204 x0 = (x2*g1 - x1*g2)/(g1 - g2);
dimensionedScalar sign(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
label vofCutCell(const label celli, const scalar alpha1, const scalar tol, const label maxIter, vector normal)
Finds matching cutValue for the given value fraction.
label calcSubCell(const label celli, const scalar cutValue, const vector &normal)
Sets internal values and returns face status.
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
const dimensionedScalar e
Elementary charge.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
constexpr scalar pi(M_PI)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
surfaceIteratorPLIC(const fvMesh &mesh, const scalar tol)
Construct from fvMesh and a scalarField.
#define WarningInFunction
Report a warning using Foam::Warning.
Mesh data needed to do the Finite Volume discretisation.
scalar VolumeOfFluid() const noexcept
Returns volume of fluid value.
const volVectorField & C() const
Return cell centres as volVectorField.
List< label > labelList
A List of labels.
SquareMatrix< scalar > scalarSquareMatrix
const labelListList & cellPoints() const
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
Solve the matrix using LU decomposition with pivoting returning the LU form of the matrix and the sol...
const volScalarField & alpha1