59 <<
"normal length is zero in cell: " << celli <<
nl
60 <<
"try increasing nCorrectors" <<
endl;
68 const labelList& pLabels = mesh_.cellPoints(celli);
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]];
101 cutCell_.calcSubCell(celli, f3, normal);
102 a3 = cutCell_.VolumeOfFluid();
105 L1 = L3; f1 = f3; a1 = a3;
109 L2 = L3; f2 = f3; a2 = a3;
113 return cutCell_.calcSubCell(celli, f3, normal);
117 if (
mag(f1 - f2) < 10*SMALL)
119 return cutCell_.calcSubCell(celli, f1, normal);
122 if (
mag(a1 - a2) < tol)
124 return cutCell_.calcSubCell(celli, 0.5*(f1 + f2), normal);
133 f3 = f1 + (f2 - f1)/3;
134 cutCell_.calcSubCell(celli, f3, normal);
135 a3 = cutCell_.VolumeOfFluid();
137 scalar f4 = f1 + 2*(f2 - f1)/3;
138 cutCell_.calcSubCell(celli, f4, normal);
139 scalar a4 = cutCell_.VolumeOfFluid();
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);
178 const scalar VOF = cutCell_.VolumeOfFluid();
194 scalar x1 =
max(1
e-3*(f2 - f1), 100*SMALL);
197 cutCell_.calcSubCell(celli, x1,normal);
198 scalar g1 = cutCell_.VolumeOfFluid() -
alpha1;
202 while (res > tol && nIter < maxIter && g1 != g2)
204 x0 = (x2*g1 - x1*g2)/(g1 - g2);
205 status = cutCell_.calcSubCell(celli, x0, normal);
206 g0 = cutCell_.VolumeOfFluid() -
alpha1;
const volScalarField & alpha1
Graphite solid properties.
T & first()
Return the first element of the list.
void size(const label n)
Older name for setAddressableSize.
T & last()
Return the last element of the list.
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
Mesh data needed to do the Finite Volume discretisation.
Finds the cutValue that matches the volume fraction.
label vofCutCell(const label celli, const scalar alpha1, const scalar tol, const label maxIter, vector normal)
Finds matching cutValue for the given value fraction.
#define WarningInFunction
Report a warning using Foam::Warning.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedScalar sign(const dimensionedScalar &ds)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
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 min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.