49 for (
const label facei : faceIDs)
66 solveScalar sumA = 0.0;
84 solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint);
85 solveScalar a =
mag(
n);
93 if (sumA < ROOTVSMALL)
95 fCtrs[facei] = fCentre;
100 fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
101 fAreas[facei] = 0.5*sumN;
110 const primitiveMesh&
mesh,
114 const List<cell>&
cells,
121 PrecisionAdaptor<solveVector, vector> tcellCtrs(cellCtrs_s,
false);
122 PrecisionAdaptor<solveScalar, scalar> tcellVols(cellVols_s,
false);
123 Field<solveVector>& cellCtrs = tcellCtrs.ref();
124 Field<solveScalar>&
cellVols = tcellVols.ref();
142 const label celli =
cellIDs[i];
143 const cell&
c =
cells[celli];
144 for (
const auto facei :
c)
146 Cc0[i] += fCtrs[facei];
153 Cc0[i] /= nCellFaces[i];
158 for (
const label celli :
cellIDs)
160 cellCtrs[celli] =
Zero;
168 const label celli =
cellIDs[i];
169 const auto&
c =
cells[celli];
172 for (
const label facei :
c)
177 const solveScalar pyr3Vol = own[facei] == celli
182 const solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cc;
185 cellCtrs[celli] += pyr3Vol*pc;
198 cellCtrs[celli] = Cc0[i];
206 const UList<face>& fcs,
213 fCtrs.resize_nocopy(fcs.size());
214 fAreas.resize_nocopy(fcs.size());
218 const face&
f = fcs[facei];
231 solveScalar sumA =
Zero;
247 solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint);
248 solveScalar a =
mag(
n);
257 if (sumA < ROOTVSMALL)
259 fCtrs[facei] = fCentre;
260 fAreas[facei] =
Zero;
264 fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
265 fAreas[facei] = 0.5*sumN;
274 const primitiveMesh&
mesh,
280 makeFaceCentresAndAreas(
mesh.
faces(),
p, fCtrs, fAreas);
314 ++nCellFaces[own[facei]];
320 ++nCellFaces[nei[facei]];
325 cEst[celli] /= nCellFaces[celli];
334 solveScalar pyr3Vol = fA & (fc - cEst[own[facei]]);
337 solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[own[facei]];
340 cellCtrs[own[facei]] += pyr3Vol*pc;
352 solveScalar pyr3Vol = fA & (cEst[nei[facei]] - fc);
355 solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[nei[facei]];
358 cellCtrs[nei[facei]] += pyr3Vol*pc;
372 cellCtrs[celli] = cEst[celli];
382 const UList<face>& fcs,
392 const face&
f = fcs[facei];
394 vector Cpf = fCtrs[facei] - ownCc;
400 - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
401 vector svHat = sv/(
mag(sv) + ROOTVSMALL);
406 scalar fd = 0.2*
mag(d) + ROOTVSMALL;
409 fd =
max(fd,
mag(svHat & (
p[
f[
pi]] - fCtrs[facei])));
428 const face&
f = fcs[facei];
430 vector Cpf = fCtrs[facei] - ownCc;
433 vector d = normal*(normal & Cpf);
438 - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
439 vector svHat = sv/(
mag(sv) + ROOTVSMALL);
444 scalar fd = 0.4*
mag(d) + ROOTVSMALL;
447 fd =
max(fd,
mag(svHat & (
p[
f[
pi]] - fCtrs[facei])));
467 return faceSkewness(
mesh.
faces(),
p, fCtrs, fAreas, facei, ownCc, neiCc);
482 return boundaryFaceSkewness(
mesh.
faces(),
p, fCtrs, fAreas, facei, ownCc);
512 auto& ortho = tortho.ref();
517 ortho[facei] = faceOrthogonality
531 const primitiveMesh&
mesh,
543 auto&
skew = tskew.ref();
548 skew[facei] = faceSkewness
556 cellCtrs[own[facei]],
567 skew[facei] = boundaryFaceSkewness
584 const primitiveMesh&
mesh,
623 const primitiveMesh&
mesh,
624 const Vector<label>& meshD,
644 sumClosed[own[facei]] += areas[facei];
645 sumMagClosed[own[facei]] +=
cmptMag(areas[facei]);
651 sumClosed[nei[facei]] -= areas[facei];
652 sumMagClosed[nei[facei]] +=
cmptMag(areas[facei]);
672 scalar maxOpenness = 0;
679 mag(sumClosed[celli][cmpt])
680 /(sumMagClosed[celli][cmpt] + ROOTVSMALL)
683 openness[celli] = maxOpenness;
687 scalar minCmpt = VGREAT;
688 scalar maxCmpt = -VGREAT;
693 minCmpt =
min(minCmpt, sumMagClosed[celli][dir]);
694 maxCmpt =
max(maxCmpt, sumMagClosed[celli][dir]);
698 scalar aspectRatio = maxCmpt/(minCmpt + ROOTVSMALL);
701 scalar v =
max(ROOTVSMALL, vols[celli]);
706 1.0/6.0*
cmptSum(sumMagClosed[celli])/
pow(v, 2.0/3.0)
710 aratio[celli] = aspectRatio;
718 const primitiveMesh&
mesh,
726 faceNormals /=
mag(faceNormals) + ROOTVSMALL;
729 auto&& faceAngles = tfaceAngles.ref();
733 const face&
f = fcs[facei];
737 scalar magEPrev =
mag(ePrev);
738 ePrev /= magEPrev + ROOTVSMALL;
740 scalar maxEdgeSin = 0.0;
745 vector e10(
p[
f.nextLabel(fp0)] -
p[
f.thisLabel(fp0)]);
746 scalar magE10 =
mag(e10);
747 e10 /= magE10 + ROOTVSMALL;
749 if (magEPrev > SMALL && magE10 > SMALL)
751 vector edgeNormal = ePrev ^ e10;
752 scalar magEdgeNormal =
mag(edgeNormal);
754 if (magEdgeNormal < maxSin)
761 edgeNormal /= magEdgeNormal;
763 if ((edgeNormal & faceNormals[facei]) < SMALL)
765 maxEdgeSin =
max(maxEdgeSin, magEdgeNormal);
774 faceAngles[facei] = maxEdgeSin;
783 const primitiveMesh&
mesh,
796 auto& faceFlatness = tfaceFlatness.ref();
800 const face&
f = fcs[facei];
802 if (
f.
size() > 3 && magAreas[facei] > ROOTVSMALL)
809 solveScalar sumA = 0.0;
817 solveVector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
821 faceFlatness[facei] = magAreas[facei]/(sumA + ROOTVSMALL);
825 return tfaceFlatness;
831 const primitiveMesh&
mesh,
832 const Vector<label>& meshD,
834 const bitSet& internalOrCoupledFace
853 auto& cellDeterminant = tcellDeterminant.ref();
859 cellDeterminant = 1.0;
870 label nInternalFaces = 0;
874 if (internalOrCoupledFace.test(curFaces[i]))
876 avgArea +=
mag(faceAreas[curFaces[i]]);
882 if (nInternalFaces == 0 || avgArea < ROOTVSMALL)
884 cellDeterminant[celli] = 0;
888 avgArea /= nInternalFaces;
894 if (internalOrCoupledFace.test(curFaces[i]))
896 areaTensor +=
sqr(faceAreas[curFaces[i]]/avgArea);
926 cellDeterminant[celli] =
mag(
det(areaTensor))/8.0;
931 return tcellDeterminant;
Field< label > labelField
Specialisation of Field<T> for label.
void size(const label n)
Older name for setAddressableSize.
List< cell > cellList
List of cell.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
Return the sum of components of a SphericalTensor.
virtual const labelList & faceNeighbour() const
Return face neighbour.
dimensionedTensor skew(const dimensionedTensor &dt)
Cell-face mesh analysis engine.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
T & first()
Access first element of the list, position [0].
const cellList & cells() const
::Foam::direction nComponents(const expressions::valueTypeCode) noexcept
The number of components associated with given valueTypeCode.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
A non-const Field/List wrapper with possible data conversion.
label nFaces() const noexcept
Number of mesh faces.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
#define forAll(list, i)
Loop across all elements in list.
pyramid< point, const point &, const face & > pyramidPointFaceRef
A pyramid using referred point and face.
Vector< solveScalar > solveVector
List< face > faceList
List of faces.
vectorField pointField
pointField is a vectorField.
void setSize(const label n)
Alias for resize()
SymmTensor< scalar > symmTensor
SymmTensor of scalars, i.e. SymmTensor<scalar>.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual const labelList & faceOwner() const
Return face owner.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
Point centre() const
Return centre (centroid)
constexpr scalar pi(M_PI)
label nInternalFaces() const noexcept
Number of internal faces.
virtual const faceList & faces() const
Return raw faces.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Templated 3D Vector derived from VectorSpace adding construction from 3 components, element access using x(), y() and z() member functions and the inner-product (dot-product) and cross-product operators.
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
const scalarField & cellVols
T & last()
Access last element of the list, position [size()-1].
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
vector point
Point is a vector.
label nCells() const noexcept
Number of mesh cells.
const dimensionedScalar c
Speed of light in a vacuum.
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< label > labelList
A List of labels.
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))
vector areaNormal() const
The area normal - with magnitude equal to area of triangle.
static constexpr const zero Zero
Global zero (0)