38 namespace reconstruction
48 void Foam::reconstruction::plicRDF::interpolateNormal()
52 leastSquareGrad<scalar> lsGrad(
"polyDegree1",mesh_.
geometricD());
74 DynamicField<vector > cellCentre(100);
75 DynamicField<scalar > alphaValues(100);
77 DynamicList<vector> foundNormals(30);
89 const label gblIdx = stencil[celli][i];
99 estimatedNormal +=
n /
max(
mag(distanceToIntSeg), SMALL);
100 weight += 1/
max(
mag(distanceToIntSeg), SMALL);
101 foundNormals.append(
n);
105 if (weight != 0 &&
mag(estimatedNormal) != 0)
107 estimatedNormal /= weight;
108 estimatedNormal /=
mag(estimatedNormal);
111 bool tooCoarse =
false;
113 if (foundNormals.size() > 1 &&
mag(estimatedNormal) != 0)
119 if ((estimatedNormal & foundNormals[i]) <= 0.98)
132 if (
mag(estimatedNormal) != 0 && !tooCoarse)
134 interfaceNormal_[i] = estimatedNormal;
143 const label gblIdx = stencil[celli][i];
146 exchangeFields_.
getValue(mesh_.
C(), mapCC, gblIdx)
153 cellCentre -= mesh_.
C()[celli];
154 interfaceNormal_[i] = lsGrad.grad(cellCentre, alphaValues);
162 leastSquareGrad<scalar> lsGrad(
"polyDegree1", mesh_.geometricD());
164 exchangeFields_.setUpCommforZone(interfaceCell_,
false);
168 exchangeFields_.getDatafromOtherProc(interfaceCell_, mesh_.C())
172 exchangeFields_.getDatafromOtherProc(interfaceCell_,
phi)
175 DynamicField<vector> cellCentre(100);
176 DynamicField<scalar> phiValues(100);
180 forAll(interfaceLabels_, i)
182 const label celli = interfaceLabels_[i];
187 for (
const label gblIdx : stencil[celli])
191 exchangeFields_.getValue(mesh_.C(), mapCC, gblIdx)
195 exchangeFields_.getValue(
phi, mapPhi, gblIdx)
199 cellCentre -= mesh_.C()[celli];
200 interfaceNormal_[i] = lsGrad.grad(cellCentre, phiValues);
205 void Foam::reconstruction::plicRDF::setInitNormals(
bool interpolate)
207 interfaceLabels_.clear();
211 if (sIterPLIC_.isASurfaceCell(alpha1_[celli]))
213 interfaceCell_[celli] =
true;
214 interfaceLabels_.append(celli);
217 interfaceNormal_.setSize(interfaceLabels_.size());
219 RDF_.markCellsNearSurf(interfaceCell_, 1);
220 const boolList& nextToInterface_ = RDF_.nextToInterface();
221 exchangeFields_.updateStencil(nextToInterface_);
234 void Foam::reconstruction::plicRDF::calcResidual
236 Map<scalar>& normalResidual,
237 Map<scalar>& avgAngle
240 exchangeFields_.setUpCommforZone(interfaceCell_,
false);
242 Map<vector> mapNormal
244 exchangeFields_.getDatafromOtherProc(interfaceCell_, normal_)
249 normalResidual.
clear();
251 forAll(interfaceLabels_, i)
253 const label celli = interfaceLabels_[i];
254 if (
mag(normal_[celli]) == 0 ||
mag(interfaceNormal_[i]) == 0)
259 scalar avgDiffNormal = 0;
260 scalar maxDiffNormal = GREAT;
262 const vector cellNormal = normal_[celli]/
mag(normal_[celli]);
264 for (
const label gblIdx : stencil[celli])
267 exchangeFields_.getValue(normal_, mapNormal, gblIdx);
269 if (
mag(normal) != 0 && i != 0)
272 scalar cosAngle =
max(
min((cellNormal &
n), 1), -1);
273 avgDiffNormal +=
acos(cosAngle) *
mag(normal);
274 weight +=
mag(normal);
275 if (cosAngle < maxDiffNormal)
277 maxDiffNormal = cosAngle;
284 avgDiffNormal /= weight;
293 scalar normalRes = (1 - (cellNormal & newCellNormal));
294 avgAngle.insert(celli, avgDiffNormal);
295 normalResidual.insert(celli, normalRes);
302 Foam::reconstruction::plicRDF::plicRDF
320 interfaceNormal_(0.2*mesh_.nCells()),
322 isoFaceTol_(modelDict().getOrDefault<scalar>(
"isoFaceTol", 1
e-8)),
323 surfCellTol_(modelDict().getOrDefault<scalar>(
"surfCellTol", 1
e-8)),
324 tol_(modelDict().getOrDefault<scalar>(
"tol", 1
e-6)),
325 relTol_(modelDict().getOrDefault<scalar>(
"relTol", 0.1)),
326 iteration_(modelDict().getOrDefault<label>(
"iterations", 5)),
327 interpolateNormal_(modelDict().getOrDefault(
"interpolateNormal",
true)),
330 sIterPLIC_(mesh_,surfCellTol_)
332 setInitNormals(
false);
337 forAll(interfaceLabels_, i)
339 const label celli = interfaceLabels_[i];
340 if (
mag(interfaceNormal_[i]) == 0)
344 sIterPLIC_.vofCutCell
353 if (sIterPLIC_.cellStatus() == 0)
355 normal_[celli] = sIterPLIC_.surfaceArea();
356 centre_[celli] = sIterPLIC_.surfaceCentre();
357 if (
mag(normal_[celli]) == 0)
359 normal_[celli] =
Zero;
360 centre_[celli] =
Zero;
365 normal_[celli] =
Zero;
366 centre_[celli] =
Zero;
376 const bool uptodate = alreadyReconstructed(forceUpdate);
378 if (uptodate && !forceUpdate)
383 if (mesh_.topoChanging())
386 if (interfaceCell_.size() != mesh_.nCells())
388 interfaceCell_.resize(mesh_.nCells());
391 interfaceCell_ =
false;
394 setInitNormals(interpolateNormal_);
400 const boolList& nextToInterface_ = RDF_.nextToInterface();
404 for (
int iter=0; iter<iteration_; ++iter)
406 forAll(interfaceLabels_, i)
408 const label celli = interfaceLabels_[i];
409 if (
mag(interfaceNormal_[i]) == 0 || tooCoarse.
found(celli))
413 sIterPLIC_.vofCutCell
422 if (sIterPLIC_.cellStatus() == 0)
425 normal_[celli] = sIterPLIC_.surfaceArea();
426 centre_[celli] = sIterPLIC_.surfaceCentre();
427 if (
mag(normal_[celli]) == 0)
429 normal_[celli] =
Zero;
430 centre_[celli] =
Zero;
435 normal_[celli] =
Zero;
436 centre_[celli] =
Zero;
440 normal_.correctBoundaryConditions();
441 centre_.correctBoundaryConditions();
445 surfaceVectorField::Boundary nHatb(mesh_.Sf().boundaryField());
446 nHatb *= 1/(mesh_.magSf().boundaryField());
457 RDF_.updateContactAngle(alpha1_, U_, nHatb);
459 calcResidual(residual, avgAngle);
463 label resCounter = 0;
465 scalar avgNormRes = 0;
470 while (resIter.found())
472 if (avgAngleIter() > 0.26 && iter > 0)
474 tooCoarse.
set(resIter.key());
480 scalar discreteError = 0.01*
sqr(avgAngleIter());
481 if (discreteError != 0)
483 normRes= resIter()/
max(discreteError, tol_);
487 normRes= resIter()/tol_;
489 avgNormRes += normRes;
513 <<
"initial residual absolute = "
514 << avgRes/resCounter <<
nl
515 <<
"initial residual normalized = "
516 << avgNormRes/resCounter <<
nl;
522 (avgNormRes/resCounter < relTol_ || avgRes/resCounter < tol_)
525 || iter + 1 == iteration_
529 <<
"iterations = " << iter <<
nl
530 <<
"final residual absolute = "
531 << avgRes/resCounter <<
nl
532 <<
"final residual normalized = " << avgNormRes/resCounter
550 if (
mag(normal_[celli]) != 0)
552 vector n = normal_[celli]/
mag(normal_[celli]);
553 scalar cutValue = (centre_[celli] - mesh_.C()[celli]) & (
n);
560 alpha1_[celli] =
cutCell.VolumeOfFluid();
564 alpha1_.correctBoundaryConditions();
565 alpha1_.oldTime () = alpha1_;
566 alpha1_.oldTime().correctBoundaryConditions();