48 for (
const label facei : faceIDs)
57 fCtrs[facei] = (1.0/3.0)*(
p[
f[0]] +
p[
f[1]] +
p[
f[2]]);
58 fAreas[facei] = 0.5*((
p[
f[1]] -
p[
f[0]])^(
p[
f[2]] -
p[
f[0]]));
65 solveScalar sumA = 0.0;
69 for (label pi = 1; pi <
nPoints; ++pi)
76 for (label pi = 0; pi <
nPoints; ++pi)
78 const label nextPi(pi ==
nPoints-1 ? 0 : pi+1);
83 solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint);
84 solveScalar a =
mag(
n);
92 if (sumA < ROOTVSMALL)
94 fCtrs[facei] = fCentre;
99 fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
100 fAreas[facei] = 0.5*sumN;
141 const label celli =
cellIDs[i];
143 for (
const auto facei : c)
145 Cc0[i] += fCtrs[facei];
152 Cc0[i] /= nCellFaces[i];
157 for (
const label celli :
cellIDs)
159 cellCtrs[celli] =
Zero;
167 const label celli =
cellIDs[i];
168 const auto& c =
cells[celli];
171 for (
const label facei : c)
176 const solveScalar pyr3Vol = own[facei] == celli
181 const solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cc;
184 cellCtrs[celli] += pyr3Vol*pc;
197 cellCtrs[celli] = Cc0[i];
217 const face&
f = fcs[facei];
224 fCtrs[facei] = (1.0/3.0)*(
p[
f[0]] +
p[
f[1]] +
p[
f[2]]);
225 fAreas[facei] = 0.5*((
p[
f[1]] -
p[
f[0]])^(
p[
f[2]] -
p[
f[0]]));
230 solveScalar sumA =
Zero;
234 for (label pi = 1; pi <
nPoints; ++pi)
240 for (label pi = 0; pi <
nPoints; ++pi)
246 solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint);
247 solveScalar a =
mag(
n);
256 if (sumA < ROOTVSMALL)
258 fCtrs[facei] = fCentre;
259 fAreas[facei] =
Zero;
263 fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
264 fAreas[facei] = 0.5*sumN;
279 makeFaceCentresAndAreas(
mesh.
faces(),
p, fCtrs, fAreas);
313 ++nCellFaces[own[facei]];
319 ++nCellFaces[nei[facei]];
324 cEst[celli] /= nCellFaces[celli];
333 solveScalar pyr3Vol = fA & (fc - cEst[own[facei]]);
336 solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[own[facei]];
339 cellCtrs[own[facei]] += pyr3Vol*pc;
351 solveScalar pyr3Vol = fA & (cEst[nei[facei]] - fc);
354 solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[nei[facei]];
357 cellCtrs[nei[facei]] += pyr3Vol*pc;
371 cellCtrs[celli] = cEst[celli];
391 const face&
f = fcs[facei];
393 vector Cpf = fCtrs[facei] - ownCc;
399 - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
400 vector svHat = sv/(
mag(sv) + ROOTVSMALL);
405 scalar fd = 0.2*
mag(d) + ROOTVSMALL;
408 fd =
max(fd,
mag(svHat & (
p[
f[pi]] - fCtrs[facei])));
427 const face&
f = fcs[facei];
429 vector Cpf = fCtrs[facei] - ownCc;
432 vector d = normal*(normal & Cpf);
437 - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
438 vector svHat = sv/(
mag(sv) + ROOTVSMALL);
443 scalar fd = 0.4*
mag(d) + ROOTVSMALL;
446 fd =
max(fd,
mag(svHat & (
p[
f[pi]] - fCtrs[facei])));
466 return faceSkewness(
mesh.
faces(),
p, fCtrs, fAreas, facei, ownCc, neiCc);
481 return boundaryFaceSkewness(
mesh.
faces(),
p, fCtrs, fAreas, facei, ownCc);
494 return (d &
s)/(
mag(d)*
mag(
s) + ROOTVSMALL);
511 auto& ortho = tortho.ref();
516 ortho[facei] = faceOrthogonality
542 auto&
skew = tskew.ref();
547 skew[facei] = faceSkewness
555 cellCtrs[own[facei]],
566 skew[facei] = boundaryFaceSkewness
643 sumClosed[own[facei]] += areas[facei];
644 sumMagClosed[own[facei]] +=
cmptMag(areas[facei]);
650 sumClosed[nei[facei]] -= areas[facei];
651 sumMagClosed[nei[facei]] +=
cmptMag(areas[facei]);
671 scalar maxOpenness = 0;
678 mag(sumClosed[celli][cmpt])
679 /(sumMagClosed[celli][cmpt] + ROOTVSMALL)
682 openness[celli] = maxOpenness;
686 scalar minCmpt = VGREAT;
687 scalar maxCmpt = -VGREAT;
692 minCmpt =
min(minCmpt, sumMagClosed[celli][dir]);
693 maxCmpt =
max(maxCmpt, sumMagClosed[celli][dir]);
697 scalar aspectRatio = maxCmpt/(minCmpt + ROOTVSMALL);
700 scalar v =
max(ROOTVSMALL, vols[celli]);
705 1.0/6.0*
cmptSum(sumMagClosed[celli])/
pow(v, 2.0/3.0)
709 aratio[celli] = aspectRatio;
728 auto&& faceAngles = tfaceAngles.ref();
732 const face&
f = fcs[facei];
736 scalar magEPrev =
mag(ePrev);
737 ePrev /= magEPrev + ROOTVSMALL;
739 scalar maxEdgeSin = 0.0;
744 vector e10(
p[
f.nextLabel(fp0)] -
p[
f.thisLabel(fp0)]);
745 scalar magE10 =
mag(e10);
746 e10 /= magE10 + ROOTVSMALL;
748 if (magEPrev > SMALL && magE10 > SMALL)
750 vector edgeNormal = ePrev ^ e10;
751 scalar magEdgeNormal =
mag(edgeNormal);
753 if (magEdgeNormal < maxSin)
760 edgeNormal /= magEdgeNormal;
764 maxEdgeSin =
max(maxEdgeSin, magEdgeNormal);
773 faceAngles[facei] = maxEdgeSin;
795 auto& faceFlatness = tfaceFlatness.ref();
799 const face&
f = fcs[facei];
801 if (
f.
size() > 3 && magAreas[facei] > ROOTVSMALL)
808 solveScalar sumA = 0.0;
816 solveVector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
820 faceFlatness[facei] = magAreas[facei]/(sumA + ROOTVSMALL);
824 return tfaceFlatness;
833 const bitSet& internalOrCoupledFace
852 auto& cellDeterminant = tcellDeterminant.ref();
858 cellDeterminant = 1.0;
869 label nInternalFaces = 0;
873 if (internalOrCoupledFace.
test(curFaces[i]))
875 avgArea +=
mag(faceAreas[curFaces[i]]);
881 if (nInternalFaces == 0 || avgArea < ROOTVSMALL)
883 cellDeterminant[celli] = 0;
887 avgArea /= nInternalFaces;
893 if (internalOrCoupledFace.
test(curFaces[i]))
895 areaTensor +=
sqr(faceAreas[curFaces[i]]/avgArea);
925 cellDeterminant[celli] =
mag(
det(areaTensor))/8.0;
930 return tcellDeterminant;
void setSize(const label n)
Alias for resize()
void resize_nocopy(const label len)
Adjust allocated size of list without necessarily.
A non-const Field/List wrapper with possible data conversion.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
T & first()
Return the first element of the list.
void size(const label n)
Older name for setAddressableSize.
T & last()
Return the last element of the list.
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
tmp< scalarField > faceSkewness() const
Return face skewness.
A cell is defined as a list of faces with extra functionality.
A face is a list of labels corresponding to mesh vertices.
static constexpr direction nComponents
Number of components in bool is 1.
virtual const faceList & faces() const
Return raw faces.
virtual const labelList & faceOwner() const
Return face owner.
virtual const labelList & faceNeighbour() const
Return face neighbour.
Cell-face mesh analysis engine.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
label nInternalFaces() const noexcept
Number of internal faces.
label nCells() const noexcept
Number of mesh cells.
label nFaces() const noexcept
Number of mesh faces.
const cellList & cells() const
A class for managing temporary objects.
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))
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf())
pyramid< point, const point &, const face & > pyramidPointFaceRef
A pyramid using referred points and faces.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
Return the sum of components of a SphericalTensor.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
static constexpr const zero Zero
Global zero (0)
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
Vector< solveScalar > solveVector
dimensionedTensor skew(const dimensionedTensor &dt)
#define forAll(list, i)
Loop across all elements in list.
const scalarField & cellVols