39 namespace reconstruction
49 void Foam::reconstruction::plicRDF::interpolateNormal()
55 leastSquareGrad<scalar> lsGrad(
"polyDegree1",mesh_.
geometricD());
77 DynamicField<vector > cellCentre(100);
78 DynamicField<scalar > alphaValues(100);
80 DynamicList<vector> foundNormals(30);
90 for (
const label gblIdx : stencil[celli])
93 exchangeFields.getValue(
normal_, mapNormal, gblIdx);
99 exchangeFields.getValue(
centre_, mapCentre, gblIdx);
101 estimatedNormal +=
n /
max(
mag(distanceToIntSeg), SMALL);
102 weight += 1/
max(
mag(distanceToIntSeg), SMALL);
103 foundNormals.append(
n);
107 if (weight != 0 &&
mag(estimatedNormal) != 0)
109 estimatedNormal /= weight;
110 estimatedNormal /=
mag(estimatedNormal);
113 bool tooCoarse =
false;
115 if (foundNormals.size() > 1 &&
mag(estimatedNormal) != 0)
121 if ((estimatedNormal & foundNormals[i]) <= 0.98)
134 if (
mag(estimatedNormal) != 0 && !tooCoarse)
136 interfaceNormal_[i] = estimatedNormal;
145 const label gblIdx = stencil[celli][i];
148 exchangeFields.getValue(mesh_.
C(), mapCC, gblIdx)
152 exchangeFields.getValue(
alpha1_, mapAlpha, gblIdx)
155 cellCentre -= mesh_.
C()[celli];
156 interfaceNormal_[i] = lsGrad.grad(cellCentre, alphaValues);
165 leastSquareGrad<scalar> lsGrad(
"polyDegree1", mesh_.geometricD());
168 exchangeFields.setUpCommforZone(interfaceCell_,
false);
172 exchangeFields.getDatafromOtherProc(interfaceCell_, mesh_.C())
176 exchangeFields.getDatafromOtherProc(interfaceCell_,
phi)
179 DynamicField<vector> cellCentre(100);
180 DynamicField<scalar> phiValues(100);
184 forAll(interfaceLabels_, i)
186 const label celli = interfaceLabels_[i];
191 for (
const label gblIdx : stencil[celli])
195 exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
199 exchangeFields.getValue(
phi, mapPhi, gblIdx)
203 cellCentre -= mesh_.C()[celli];
204 interfaceNormal_[i] = lsGrad.grad(cellCentre, phiValues);
209 void Foam::reconstruction::plicRDF::setInitNormals(
bool interpolate)
214 interfaceLabels_.clear();
218 if (sIterPLIC_.isASurfaceCell(alpha1_[celli]))
220 interfaceCell_[celli] =
true;
221 interfaceLabels_.append(celli);
224 interfaceNormal_.setSize(interfaceLabels_.size());
226 RDF_.markCellsNearSurf(interfaceCell_, 1);
227 const boolList& nextToInterface_ = RDF_.nextToInterface();
228 exchangeFields.updateStencil(nextToInterface_);
241 void Foam::reconstruction::plicRDF::calcResidual
243 List<normalRes>& normalResidual
248 exchangeFields.setUpCommforZone(interfaceCell_,
false);
250 Map<vector> mapNormal
252 exchangeFields.getDatafromOtherProc(interfaceCell_, normal_)
257 forAll(interfaceLabels_, i)
259 const label celli = interfaceLabels_[i];
260 if (
mag(normal_[celli]) == 0 ||
mag(interfaceNormal_[i]) == 0)
262 normalResidual[i].celli = celli;
263 normalResidual[i].normalResidual = 0;
264 normalResidual[i].avgAngle = 0;
268 scalar avgDiffNormal = 0;
269 scalar maxDiffNormal = GREAT;
271 const vector cellNormal = normal_[celli]/
mag(normal_[celli]);
274 const label gblIdx = stencil[celli][j];
276 exchangeFields.getValue(normal_, mapNormal, gblIdx);
278 if (
mag(normal) != 0 && j != 0)
281 scalar cosAngle =
max(
min((cellNormal &
n), 1), -1);
282 avgDiffNormal +=
acos(cosAngle) *
mag(normal);
283 weight +=
mag(normal);
284 if (cosAngle < maxDiffNormal)
286 maxDiffNormal = cosAngle;
293 avgDiffNormal /= weight;
302 scalar normalRes = (1 - (cellNormal & newCellNormal));
303 normalResidual[i].celli = celli;
304 normalResidual[i].normalResidual = normalRes;
305 normalResidual[i].avgAngle = avgDiffNormal;
312 Foam::reconstruction::plicRDF::plicRDF
330 interfaceNormal_(0.2*mesh_.nCells()),
332 isoFaceTol_(modelDict().getOrDefault<scalar>(
"isoFaceTol", 1
e-8)),
333 surfCellTol_(modelDict().getOrDefault<scalar>(
"surfCellTol", 1
e-8)),
334 tol_(modelDict().getOrDefault<scalar>(
"tol", 1
e-6)),
335 relTol_(modelDict().getOrDefault<scalar>(
"relTol", 0.1)),
336 iteration_(modelDict().getOrDefault<label>(
"iterations", 5)),
337 interpolateNormal_(modelDict().getOrDefault(
"interpolateNormal",
true)),
339 sIterPLIC_(mesh_,surfCellTol_)
341 setInitNormals(
false);
346 forAll(interfaceLabels_, i)
348 const label celli = interfaceLabels_[i];
349 if (
mag(interfaceNormal_[i]) == 0)
353 sIterPLIC_.vofCutCell
362 if (sIterPLIC_.cellStatus() == 0)
364 normal_[celli] = sIterPLIC_.surfaceArea();
365 centre_[celli] = sIterPLIC_.surfaceCentre();
366 if (
mag(normal_[celli]) == 0)
368 normal_[celli] =
Zero;
369 centre_[celli] =
Zero;
374 normal_[celli] =
Zero;
375 centre_[celli] =
Zero;
387 const bool uptodate = alreadyReconstructed(forceUpdate);
389 if (uptodate && !forceUpdate)
394 if (mesh_.topoChanging())
397 if (interfaceCell_.size() != mesh_.nCells())
399 interfaceCell_.resize(mesh_.nCells());
402 interfaceCell_ =
false;
405 setInitNormals(interpolateNormal_);
411 const boolList& nextToInterface_ = RDF_.nextToInterface();
413 bitSet tooCoarse(mesh_.nCells(),
false);
415 for (label iter=0; iter<iteration_; ++iter)
417 forAll(interfaceLabels_, i)
419 const label celli = interfaceLabels_[i];
420 if (
mag(interfaceNormal_[i]) == 0 || tooCoarse.test(celli))
424 sIterPLIC_.vofCutCell
433 if (sIterPLIC_.cellStatus() == 0)
436 normal_[celli] = sIterPLIC_.surfaceArea();
437 centre_[celli] = sIterPLIC_.surfaceCentre();
438 if (
mag(normal_[celli]) == 0)
440 normal_[celli] =
Zero;
441 centre_[celli] =
Zero;
446 normal_[celli] =
Zero;
447 centre_[celli] =
Zero;
451 normal_.correctBoundaryConditions();
452 centre_.correctBoundaryConditions();
455 surfaceVectorField::Boundary nHatb(mesh_.Sf().boundaryField());
456 nHatb *= 1/(mesh_.magSf().boundaryField());
467 RDF_.updateContactAngle(alpha1_, U_, nHatb);
469 calcResidual(normalResidual);
472 label resCounter = 0;
474 scalar avgNormRes = 0;
479 const label celli = normalResidual[i].celli;
480 const scalar normalRes= normalResidual[i].normalResidual;
481 const scalar avgA = normalResidual[i].avgAngle;
483 if (avgA > 0.26 && iter > 0)
485 tooCoarse.set(celli);
491 scalar discreteError = 0.01*
sqr(avgA);
492 if (discreteError != 0)
494 normRes= normalRes/
max(discreteError, tol_);
498 normRes= normalRes/tol_;
500 avgNormRes += normRes;
521 <<
"initial residual absolute = "
522 << avgRes/resCounter <<
nl
523 <<
"initial residual normalized = "
524 << avgNormRes/resCounter <<
nl;
530 (avgNormRes/resCounter < relTol_ || avgRes/resCounter < tol_)
533 || iter + 1 == iteration_
537 <<
"iterations = " << iter <<
nl
538 <<
"final residual absolute = "
539 << avgRes/resCounter <<
nl
540 <<
"final residual normalized = " << avgNormRes/resCounter
558 if (
mag(normal_[celli]) != 0)
560 vector n = normal_[celli]/
mag(normal_[celli]);
561 scalar cutValue = (centre_[celli] - mesh_.C()[celli]) & (
n);
568 alpha1_[celli] =
cutCell.VolumeOfFluid();
572 alpha1_.correctBoundaryConditions();
573 alpha1_.oldTime () = alpha1_;
574 alpha1_.oldTime().correctBoundaryConditions();