40 namespace reconstruction
50 void Foam::reconstruction::plicRDF::interpolateNormal()
55 leastSquareGrad<scalar> lsGrad(
"polyDegree1",mesh_.
geometricD());
78 DynamicField<vector> cellCentre(100);
79 DynamicField<scalar> alphaValues(100);
81 DynamicList<vector> foundNormals(30);
91 for (
const label gblIdx : stencil[celli])
94 exchangeFields.getValue(
normal_, mapNormal, gblIdx);
100 exchangeFields.getValue(
centre_, mapCentre, gblIdx);
102 estimatedNormal +=
n /
max(
mag(distanceToIntSeg), SMALL);
103 weight += 1/
max(
mag(distanceToIntSeg), SMALL);
104 foundNormals.append(
n);
108 if (weight != 0 &&
mag(estimatedNormal) != 0)
110 estimatedNormal /= weight;
111 estimatedNormal /=
mag(estimatedNormal);
114 bool tooCoarse =
false;
116 if (foundNormals.size() > 1 &&
mag(estimatedNormal) != 0)
122 if ((estimatedNormal & foundNormals[i]) <= 0.98)
135 if (
mag(estimatedNormal) != 0 && !tooCoarse)
137 interfaceNormal_[i] = estimatedNormal;
146 const label gblIdx = stencil[celli][i];
149 exchangeFields.getValue(mesh_.
C(), mapCC, gblIdx)
153 exchangeFields.getValue(
alpha1_, mapAlpha, gblIdx)
156 cellCentre -= mesh_.
C()[celli];
157 interfaceNormal_[i] = lsGrad.grad(cellCentre, alphaValues);
166 leastSquareGrad<scalar> lsGrad(
"polyDegree1", mesh_.geometricD());
169 exchangeFields.setUpCommforZone(interfaceCell_,
false);
173 exchangeFields.getDatafromOtherProc(interfaceCell_, mesh_.C())
177 exchangeFields.getDatafromOtherProc(interfaceCell_,
phi)
180 DynamicField<vector> cellCentre(100);
181 DynamicField<scalar> phiValues(100);
185 forAll(interfaceLabels_, i)
187 const label celli = interfaceLabels_[i];
192 for (
const label gblIdx : stencil[celli])
196 exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
200 exchangeFields.getValue(
phi, mapPhi, gblIdx)
204 cellCentre -= mesh_.C()[celli];
205 interfaceNormal_[i] = lsGrad.grad(cellCentre, phiValues);
210 void Foam::reconstruction::plicRDF::setInitNormals(
bool interpolate)
215 interfaceLabels_.clear();
219 if (sIterPLIC_.isASurfaceCell(alpha1_[celli]))
221 interfaceCell_[celli] =
true;
222 interfaceLabels_.append(celli);
225 interfaceNormal_.setSize(interfaceLabels_.size());
227 RDF_.markCellsNearSurf(interfaceCell_, 1);
228 const boolList& nextToInterface_ = RDF_.nextToInterface();
229 exchangeFields.updateStencil(nextToInterface_);
242 void Foam::reconstruction::plicRDF::calcResidual
244 List<normalRes>& normalResidual
249 exchangeFields.setUpCommforZone(interfaceCell_,
false);
251 Map<vector> mapNormal
253 exchangeFields.getDatafromOtherProc(interfaceCell_, normal_)
258 forAll(interfaceLabels_, i)
260 const label celli = interfaceLabels_[i];
261 if (
mag(normal_[celli]) == 0 ||
mag(interfaceNormal_[i]) == 0)
263 normalResidual[i].celli = celli;
264 normalResidual[i].normalResidual = 0;
265 normalResidual[i].avgAngle = 0;
269 scalar avgDiffNormal = 0;
270 scalar maxDiffNormal = GREAT;
272 const vector cellNormal = normal_[celli]/
mag(normal_[celli]);
275 const label gblIdx = stencil[celli][j];
277 exchangeFields.getValue(normal_, mapNormal, gblIdx);
279 if (
mag(normal) != 0 && j != 0)
282 scalar cosAngle =
max(
min((cellNormal &
n), 1), -1);
283 avgDiffNormal +=
acos(cosAngle) *
mag(normal);
284 weight +=
mag(normal);
285 if (cosAngle < maxDiffNormal)
287 maxDiffNormal = cosAngle;
294 avgDiffNormal /= weight;
303 scalar normalRes = (1 - (cellNormal & newCellNormal));
304 normalResidual[i].celli = celli;
305 normalResidual[i].normalResidual = normalRes;
306 normalResidual[i].avgAngle = avgDiffNormal;
313 Foam::reconstruction::plicRDF::plicRDF
331 interfaceNormal_(0.2*mesh_.nCells()),
333 isoFaceTol_(modelDict().getOrDefault<scalar>(
"isoFaceTol", 1
e-8)),
334 surfCellTol_(modelDict().getOrDefault<scalar>(
"surfCellTol", 1
e-8)),
335 tol_(modelDict().getOrDefault<scalar>(
"tol", 1
e-6)),
336 relTol_(modelDict().getOrDefault<scalar>(
"relTol", 0.1)),
337 iteration_(modelDict().getOrDefault<label>(
"iterations", 5)),
338 interpolateNormal_(modelDict().getOrDefault(
"interpolateNormal", true)),
340 sIterPLIC_(mesh_,surfCellTol_)
342 setInitNormals(
false);
350 if (
mag(interfaceNormal_[i]) == 0)
388 const bool uptodate = alreadyReconstructed(forceUpdate);
390 if (uptodate && !forceUpdate)
395 if (mesh_.topoChanging())
398 if (interfaceCell_.size() != mesh_.nCells())
400 interfaceCell_.resize(mesh_.nCells());
403 interfaceCell_ =
false;
406 setInitNormals(interpolateNormal_);
412 const boolList& nextToInterface_ = RDF_.nextToInterface();
414 bitSet tooCoarse(mesh_.nCells(),
false);
416 for (label iter=0; iter<iteration_; ++iter)
418 forAll(interfaceLabels_, i)
420 const label celli = interfaceLabels_[i];
421 if (
mag(interfaceNormal_[i]) == 0 || tooCoarse.test(celli))
425 sIterPLIC_.vofCutCell
434 if (sIterPLIC_.cellStatus() == 0)
437 normal_[celli] = sIterPLIC_.surfaceArea();
438 centre_[celli] = sIterPLIC_.surfaceCentre();
439 if (
mag(normal_[celli]) == 0)
441 normal_[celli] =
Zero;
442 centre_[celli] =
Zero;
447 normal_[celli] =
Zero;
448 centre_[celli] =
Zero;
452 normal_.correctBoundaryConditions();
453 centre_.correctBoundaryConditions();
454 List<normalRes> normalResidual(interfaceLabels_.size());
457 nHatb *= 1/(mesh_.magSf().boundaryField());
468 RDF_.updateContactAngle(alpha1_, U_, nHatb);
470 calcResidual(normalResidual);
473 label resCounter = 0;
475 scalar avgNormRes = 0;
480 const label celli = normalResidual[i].celli;
481 const scalar normalRes= normalResidual[i].normalResidual;
482 const scalar avgA = normalResidual[i].avgAngle;
484 if (avgA > 0.26 && iter > 0)
486 tooCoarse.set(celli);
492 scalar discreteError = 0.01*
sqr(avgA);
493 if (discreteError != 0)
495 normRes= normalRes/
max(discreteError, tol_);
499 normRes= normalRes/tol_;
501 avgNormRes += normRes;
507 reduce(avgRes,sumOp<scalar>());
508 reduce(avgNormRes,sumOp<scalar>());
509 reduce(resCounter,sumOp<label>());
522 <<
"initial residual absolute = " 523 << avgRes/resCounter <<
nl 524 <<
"initial residual normalized = " 525 << avgNormRes/resCounter <<
nl;
531 (avgNormRes/resCounter < relTol_ || avgRes/resCounter < tol_)
534 || iter + 1 == iteration_
538 <<
"iterations = " << iter <<
nl 539 <<
"final residual absolute = " 540 << avgRes/resCounter <<
nl 541 <<
"final residual normalized = " << avgNormRes/resCounter
555 cutCellPLIC cutCell(mesh_);
559 if (
mag(normal_[celli]) != 0)
561 vector n = normal_[celli]/
mag(normal_[celli]);
562 scalar cutValue = (centre_[celli] - mesh_.C()[celli]) & (
n);
569 alpha1_[celli] = cutCell.VolumeOfFluid();
573 alpha1_.correctBoundaryConditions();
574 alpha1_.oldTime () = alpha1_;
575 alpha1_.oldTime().correctBoundaryConditions();
Original code supplied by Henning Scheufler, DLR (2019)
addToRunTimeSelectionTable(reconstructionSchemes, isoAlpha, components)
List< labelList > labelListList
A List of labelList.
dimensionedScalar acos(const dimensionedScalar &ds)
scalar deltaTValue() const noexcept
Return time step value.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
A list of keyword definitions, which are a keyword followed by a number of values (eg...
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline '\n' character (0x0a)
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
volVectorField centre_
Interface centres.
label vofCutCell(const label celli, const scalar alpha1, const scalar tol, const label maxIter, vector normal)
Finds matching cutValue for the given value fraction.
GeometricBoundaryField< vector, fvsPatchField, surfaceMesh > Boundary
Type of boundary fields.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
const Time & time() const
Return the top-level database.
const vector & surfaceArea() const
The area vector of cutting isosurface.
const point & surfaceCentre() const
The centre of cutting isosurface.
Macros for easy insertion into run-time selection tables.
boolList interfaceCell_
Is interface cell?
#define forAll(list, i)
Loop across all elements in list.
virtual void reconstruct(bool forceUpdate=true)
Reconstruct interface.
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
volVectorField normal_
Interface area normals.
static zoneDistribute & New(const fvMesh &)
Selector.
const dimensionedScalar e
Elementary charge.
static const Identity< scalar > I
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
label cellStatus()
The cellStatus.
#define DebugInfo
Report an information message using Foam::Info.
volScalarField & alpha1_
Reference to the VoF Field.
#define addProfilingInFunction(name)
Define profiling trigger with specified name and description corresponding to the compiler-defined fu...
vector point
Point is a vector.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
virtual void mapAlphaField() const
Map VoF Field in case of refinement.
defineTypeNameAndDebug(isoAlpha, 0)
DynamicField< label > interfaceLabels_
List of cell labels that have an interface.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
const volVectorField & C() const
Return cell centres as volVectorField.
const volVectorField & U_
Reference to the velocity field.
List< bool > boolList
A List of bools.
const dimensionSet dimArea(sqr(dimLength))
const volVectorField & centre() const noexcept
const-Reference to interface centres
static constexpr const zero Zero
Global zero (0)
const volScalarField & alpha1