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;