57 const labelList& pLabels = mesh_.cellPoints()[celli];
61 fvert[pi] = ap_[pLabels[pi]];
65 scalar f1 = fvert[order.
first()];
66 scalar f2 = fvert[order.
last()];
80 label L2 = fvert.
size() - 1;
87 L3 = round(0.5*(L1 + L2));
88 f3 = fvert[order[L3]];
89 cutCell_.calcSubCell(celli, f3);
90 a3 = cutCell_.VolumeOfFluid();
93 L1 = L3; f1 = f3; a1 = a3;
97 L2 = L3; f2 = f3; a2 = a3;
101 return cutCell_.calcSubCell(celli, f3);
105 if (
mag(f1 - f2) < 10*SMALL)
107 return cutCell_.calcSubCell(celli, f1);
110 if (
mag(a1 - a2) < tol)
112 return cutCell_.calcSubCell(celli, 0.5*(f1 + f2));
119 f3 = f1 + (f2 - f1)/3;
120 cutCell_.calcSubCell(celli, f3);
121 a3 = cutCell_.VolumeOfFluid();
123 scalar f4 = f1 + 2*(f2 - f1)/3;
124 cutCell_.calcSubCell(celli, f4);
125 scalar a4 = cutCell_.VolumeOfFluid();
130 a[0] = a1, a[1] = a3, a[2] = a4, a[3] = a2;
131 fOld[0] = f1, fOld[1] = f3, fOld[2] = f4, fOld[3] = f2;
132 f[0] = 0,
f[1] = (f3-f1)/(f2-f1),
f[2] = (f4-f1)/(f2-f1),
f[3] = 1;
138 M[i][j] =
pow(
f[i], 3 - j);
148 f3 =
f[1]; a3 = a[1];
151 while (res > tol && nIter < 10*maxIter)
154 /(3*
C[0]*
sqr(f3) + 2*
C[1]*f3 +
C[2]);
155 a3 =
C[0]*
pow3(f3) +
C[1]*
sqr(f3) +
C[2]*f3 +
C[3];
160 f3 = f3*(f2 - f1) + f1;
163 label status = cutCell_.calcSubCell(celli, f3);
164 const scalar VOF = cutCell_.VolumeOfFluid();
180 scalar x1 =
max(1
e-3*(f2 - f1), 100*SMALL);
183 cutCell_.calcSubCell(celli, x1);
184 scalar g1 = cutCell_.VolumeOfFluid() -
alpha1;
188 while (res > tol && nIter < maxIter && g1 != g2)
190 x0 = (x2*g1 - x1*g2)/(g1 - g2);
191 status = cutCell_.calcSubCell(celli, x0);
192 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.
Mesh data needed to do the Finite Volume discretisation.
Finds the isovalue that matches the volume fraction.
label vofCutCell(const label celli, const scalar alpha1, const scalar tol, const label maxIter)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
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)
#define forAll(list, i)
Loop across all elements in list.