Go to the documentation of this file.
56 fCtrs[facei] = (1.0/3.0)*(
p[
f[0]] +
p[
f[1]] +
p[
f[2]]);
57 fAreas[facei] = 0.5*((
p[
f[1]] -
p[
f[0]])^(
p[
f[2]] -
p[
f[0]]));
63 solveVector sumN =
Zero;
64 solveScalar sumA = 0.0;
65 solveVector sumAc =
Zero;
67 solveVector fCentre =
p[
f[0]];
70 fCentre += solveVector(
p[
f[
pi]]);
78 const solveVector nextPoint(
p[
f[nextPi]]);
79 const solveVector thisPoint(
p[
f[
pi]]);
81 solveVector
c = thisPoint + nextPoint + fCentre;
82 solveVector
n = (nextPoint - thisPoint)^(fCentre - thisPoint);
83 solveScalar a =
mag(
n);
91 if (sumA < ROOTVSMALL)
93 fCtrs[facei] = fCentre;
98 fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
99 fAreas[facei] = 0.5*sumN;
137 cEst[own[facei]] += solveVector(fCtrs[facei]);
138 ++nCellFaces[own[facei]];
143 cEst[nei[facei]] += solveVector(fCtrs[facei]);
144 ++nCellFaces[nei[facei]];
149 cEst[celli] /= nCellFaces[celli];
154 const solveVector fc(fCtrs[facei]);
155 const solveVector fA(fAreas[facei]);
158 solveScalar pyr3Vol =
159 fA & (fc - cEst[own[facei]]);
162 solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[own[facei]];
165 cellCtrs[own[facei]] += pyr3Vol*pc;
173 const solveVector fc(fCtrs[facei]);
174 const solveVector fA(fAreas[facei]);
177 solveScalar pyr3Vol =
178 fA & (cEst[nei[facei]] - fc);
181 solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[nei[facei]];
184 cellCtrs[nei[facei]] += pyr3Vol*pc;
198 cellCtrs[celli] = cEst[celli];
218 vector Cpf = fCtrs[facei] - ownCc;
224 - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
225 vector svHat = sv/(
mag(sv) + ROOTVSMALL);
230 scalar fd = 0.2*
mag(d) + ROOTVSMALL;
234 fd =
max(fd,
mag(svHat & (
p[
f[
pi]] - fCtrs[facei])));
253 vector Cpf = fCtrs[facei] - ownCc;
256 vector d = normal*(normal & Cpf);
262 - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
263 vector svHat = sv/(
mag(sv) + ROOTVSMALL);
268 scalar fd = 0.4*
mag(d) + ROOTVSMALL;
272 fd =
max(fd,
mag(svHat & (
p[
f[
pi]] - fCtrs[facei])));
289 return (d &
s)/(
mag(d)*
mag(
s) + ROOTVSMALL);
311 ortho[facei] = faceOrthogonality
340 skew[facei] = faceSkewness
348 cellCtrs[own[facei]],
357 for (label facei =
mesh.nInternalFaces(); facei <
mesh.nFaces(); facei++)
359 skew[facei] = boundaryFaceSkewness
388 ownPyrVol.setSize(
mesh.nFaces());
389 neiPyrVol.setSize(
mesh.nInternalFaces());
400 if (
mesh.isInternalFace(facei))
436 sumClosed[own[facei]] += areas[facei];
437 sumMagClosed[own[facei]] +=
cmptMag(areas[facei]);
443 sumClosed[nei[facei]] -= areas[facei];
444 sumMagClosed[nei[facei]] +=
cmptMag(areas[facei]);
449 for (
direction dir = 0; dir < vector::nComponents; dir++)
459 openness.setSize(
mesh.nCells());
460 aratio.setSize(
mesh.nCells());
464 scalar maxOpenness = 0;
466 for (
direction cmpt=0; cmpt<vector::nComponents; cmpt++)
471 mag(sumClosed[celli][cmpt])
472 /(sumMagClosed[celli][cmpt] + ROOTVSMALL)
475 openness[celli] = maxOpenness;
479 scalar minCmpt = VGREAT;
480 scalar maxCmpt = -VGREAT;
481 for (
direction dir = 0; dir < vector::nComponents; dir++)
485 minCmpt =
min(minCmpt, sumMagClosed[celli][dir]);
486 maxCmpt =
max(maxCmpt, sumMagClosed[celli][dir]);
490 scalar aspectRatio = maxCmpt/(minCmpt + ROOTVSMALL);
493 scalar v =
max(ROOTVSMALL, vols[celli]);
498 1.0/6.0*
cmptSum(sumMagClosed[celli])/
pow(v, 2.0/3.0)
502 aratio[celli] = aspectRatio;
526 const face&
f = fcs[facei];
530 scalar magEPrev =
mag(ePrev);
531 ePrev /= magEPrev + ROOTVSMALL;
533 scalar maxEdgeSin = 0.0;
538 label fp1 =
f.fcIndex(fp0);
542 scalar magE10 =
mag(e10);
543 e10 /= magE10 + ROOTVSMALL;
545 if (magEPrev > SMALL && magE10 > SMALL)
547 vector edgeNormal = ePrev ^ e10;
548 scalar magEdgeNormal =
mag(edgeNormal);
550 if (magEdgeNormal < maxSin)
557 edgeNormal /= magEdgeNormal;
561 maxEdgeSin =
max(maxEdgeSin, magEdgeNormal);
570 faceAngles[facei] = maxEdgeSin;
598 const face&
f = fcs[facei];
600 if (
f.size() > 3 && magAreas[facei] > ROOTVSMALL)
602 const solveVector fc = fCtrs[facei];
607 solveScalar sumA = 0.0;
611 const solveVector thisPoint =
p[
f[fp]];
612 const solveVector nextPoint =
p[
f.nextLabel(fp)];
615 solveVector
n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
619 faceFlatness[facei] = magAreas[facei]/(sumA + ROOTVSMALL);
623 return tfaceFlatness;
632 const bitSet& internalOrCoupledFace
638 for (
direction dir = 0; dir < vector::nComponents; dir++)
657 cellDeterminant = 1.0;
668 label nInternalFaces = 0;
672 if (internalOrCoupledFace.
test(curFaces[i]))
674 avgArea +=
mag(faceAreas[curFaces[i]]);
680 if (nInternalFaces == 0 || avgArea < ROOTVSMALL)
682 cellDeterminant[celli] = 0;
686 avgArea /= nInternalFaces;
692 if (internalOrCoupledFace.
test(curFaces[i]))
694 areaTensor +=
sqr(faceAreas[curFaces[i]]/avgArea);
724 cellDeterminant[celli] =
mag(
det(areaTensor))/8.0;
729 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())
const scalarField & cellVols
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)
A non-const Field/List wrapper with possible data conversion.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr scalar pi(M_PI)
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
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)
Container< Type > & ref() const
pyramid< point, const point &, const face & > pyramidPointFaceRef
A pyramid using referred points and faces.
Cell-face mesh analysis engine.