Go to the documentation of this file.
47 vector Cpf = fCtrs[facei] - ownCc;
53 - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
59 scalar fd = 0.2*
mag(d) + ROOTVSMALL;
63 fd =
max(fd,
mag(svHat & (
p[
f[
pi]] - fCtrs[facei])));
82 vector Cpf = fCtrs[facei] - ownCc;
85 vector d = normal*(normal & Cpf);
91 - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
97 scalar fd = 0.4*
mag(d) + ROOTVSMALL;
101 fd =
max(fd,
mag(svHat & (
p[
f[
pi]] - fCtrs[facei])));
118 return (d &
s)/(
mag(d)*
mag(
s) + ROOTVSMALL);
140 ortho[facei] = faceOrthogonality
169 skew[facei] = faceSkewness
177 cellCtrs[own[facei]],
186 for (label facei =
mesh.nInternalFaces(); facei <
mesh.nFaces(); facei++)
188 skew[facei] = boundaryFaceSkewness
217 ownPyrVol.setSize(
mesh.nFaces());
218 neiPyrVol.setSize(
mesh.nInternalFaces());
229 if (
mesh.isInternalFace(facei))
265 sumClosed[own[facei]] += areas[facei];
266 sumMagClosed[own[facei]] +=
cmptMag(areas[facei]);
272 sumClosed[nei[facei]] -= areas[facei];
273 sumMagClosed[nei[facei]] +=
cmptMag(areas[facei]);
278 for (
direction dir = 0; dir < vector::nComponents; dir++)
288 openness.setSize(
mesh.nCells());
289 aratio.setSize(
mesh.nCells());
293 scalar maxOpenness = 0;
295 for (
direction cmpt=0; cmpt<vector::nComponents; cmpt++)
300 mag(sumClosed[celli][cmpt])
301 /(sumMagClosed[celli][cmpt] + ROOTVSMALL)
304 openness[celli] = maxOpenness;
308 scalar minCmpt = VGREAT;
309 scalar maxCmpt = -VGREAT;
310 for (
direction dir = 0; dir < vector::nComponents; dir++)
314 minCmpt =
min(minCmpt, sumMagClosed[celli][dir]);
315 maxCmpt =
max(maxCmpt, sumMagClosed[celli][dir]);
319 scalar aspectRatio = maxCmpt/(minCmpt + ROOTVSMALL);
322 scalar v =
max(ROOTVSMALL, vols[celli]);
327 1.0/6.0*
cmptSum(sumMagClosed[celli])/
pow(v, 2.0/3.0)
331 aratio[celli] = aspectRatio;
355 const face&
f = fcs[facei];
359 scalar magEPrev =
mag(ePrev);
360 ePrev /= magEPrev + ROOTVSMALL;
362 scalar maxEdgeSin = 0.0;
367 label fp1 =
f.fcIndex(fp0);
371 scalar magE10 =
mag(e10);
372 e10 /= magE10 + ROOTVSMALL;
374 if (magEPrev > SMALL && magE10 > SMALL)
376 vector edgeNormal = ePrev ^ e10;
377 scalar magEdgeNormal =
mag(edgeNormal);
379 if (magEdgeNormal < maxSin)
386 edgeNormal /= magEdgeNormal;
390 maxEdgeSin =
max(maxEdgeSin, magEdgeNormal);
399 faceAngles[facei] = maxEdgeSin;
427 const face&
f = fcs[facei];
429 if (
f.size() > 3 && magAreas[facei] > ROOTVSMALL)
431 const solveVector fc = fCtrs[facei];
436 solveScalar sumA = 0.0;
440 const solveVector thisPoint =
p[
f[fp]];
441 const solveVector nextPoint =
p[
f.nextLabel(fp)];
444 solveVector
n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
448 faceFlatness[facei] = magAreas[facei]/(sumA + ROOTVSMALL);
452 return tfaceFlatness;
461 const bitSet& internalOrCoupledFace
467 for (
direction dir = 0; dir < vector::nComponents; dir++)
486 cellDeterminant = 1.0;
497 label nInternalFaces = 0;
501 if (internalOrCoupledFace.
test(curFaces[i]))
503 avgArea +=
mag(faceAreas[curFaces[i]]);
509 if (nInternalFaces == 0 || avgArea < ROOTVSMALL)
511 cellDeterminant[celli] = 0;
515 avgArea /= nInternalFaces;
521 if (internalOrCoupledFace.
test(curFaces[i]))
523 areaTensor +=
sqr(faceAreas[curFaces[i]]/avgArea);
553 cellDeterminant[celli] =
mag(
det(areaTensor))/8.0;
558 return tcellDeterminant;
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
dimensionedTensor skew(const dimensionedTensor &dt)
A class for managing temporary objects.
static constexpr const zero Zero
Global zero (0)
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
#define forAll(list, i)
Loop across all elements in list.
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf())
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
Return the sum of components of a SphericalTensor.
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr scalar pi(M_PI)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
const dimensionedScalar c
Speed of light in a vacuum.
A face is a list of labels corresponding to mesh vertices.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
pyramid< point, const point &, const face & > pyramidPointFaceRef
Cell-face mesh analysis engine.