primitiveMeshTools.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2012-2016 OpenFOAM Foundation
9  Copyright (C) 2017-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "primitiveMeshTools.H"
30 #include "primitiveMesh.H"
31 #include "syncTools.H"
32 #include "pyramidPointFaceRef.H"
33 #include "PrecisionAdaptor.H"
34 
35 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36 
38 (
39  const primitiveMesh& mesh,
40  const pointField& p,
41  vectorField& fCtrs,
42  vectorField& fAreas
43 )
44 {
45  const faceList& fs = mesh.faces();
46 
47  forAll(fs, facei)
48  {
49  const labelList& f = fs[facei];
50  const label nPoints = f.size();
51 
52  // If the face is a triangle, do a direct calculation for efficiency
53  // and to avoid round-off error-related problems
54  if (nPoints == 3)
55  {
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]]));
58  }
59  else
60  {
61  typedef Vector<solveScalar> solveVector;
62 
63  solveVector sumN = Zero;
64  solveScalar sumA = 0.0;
65  solveVector sumAc = Zero;
66 
67  solveVector fCentre = p[f[0]];
68  for (label pi = 1; pi < nPoints; pi++)
69  {
70  fCentre += solveVector(p[f[pi]]);
71  }
72 
73  fCentre /= nPoints;
74 
75  for (label pi = 0; pi < nPoints; pi++)
76  {
77  const label nextPi(pi == nPoints-1 ? 0 : pi+1);
78  const solveVector nextPoint(p[f[nextPi]]);
79  const solveVector thisPoint(p[f[pi]]);
80 
81  solveVector c = thisPoint + nextPoint + fCentre;
82  solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint);
83  solveScalar a = mag(n);
84  sumN += n;
85  sumA += a;
86  sumAc += a*c;
87  }
88 
89  // This is to deal with zero-area faces. Mark very small faces
90  // to be detected in e.g., processorPolyPatch.
91  if (sumA < ROOTVSMALL)
92  {
93  fCtrs[facei] = fCentre;
94  fAreas[facei] = Zero;
95  }
96  else
97  {
98  fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
99  fAreas[facei] = 0.5*sumN;
100  }
101  }
102  }
103 }
104 
105 
107 (
108  const primitiveMesh& mesh,
109  const vectorField& fCtrs,
110  const vectorField& fAreas,
111  vectorField& cellCtrs_s,
112  scalarField& cellVols_s
113 )
114 {
115  typedef Vector<solveScalar> solveVector;
116 
117  PrecisionAdaptor<solveVector, vector> tcellCtrs(cellCtrs_s, false);
118  PrecisionAdaptor<solveScalar, scalar> tcellVols(cellVols_s, false);
119  Field<solveVector>& cellCtrs = tcellCtrs.ref();
120  Field<solveScalar>& cellVols = tcellVols.ref();
121 
122  // Clear the fields for accumulation
123  cellCtrs = Zero;
124  cellVols = 0.0;
125 
126  const labelList& own = mesh.faceOwner();
127  const labelList& nei = mesh.faceNeighbour();
128 
129  // first estimate the approximate cell centre as the average of
130  // face centres
131 
132  Field<solveVector> cEst(mesh.nCells(), Zero);
133  labelField nCellFaces(mesh.nCells(), Zero);
134 
135  forAll(own, facei)
136  {
137  cEst[own[facei]] += solveVector(fCtrs[facei]);
138  ++nCellFaces[own[facei]];
139  }
140 
141  forAll(nei, facei)
142  {
143  cEst[nei[facei]] += solveVector(fCtrs[facei]);
144  ++nCellFaces[nei[facei]];
145  }
146 
147  forAll(cEst, celli)
148  {
149  cEst[celli] /= nCellFaces[celli];
150  }
151 
152  forAll(own, facei)
153  {
154  const solveVector fc(fCtrs[facei]);
155  const solveVector fA(fAreas[facei]);
156 
157  // Calculate 3*face-pyramid volume
158  solveScalar pyr3Vol =
159  fA & (fc - cEst[own[facei]]);
160 
161  // Calculate face-pyramid centre
162  solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[own[facei]];
163 
164  // Accumulate volume-weighted face-pyramid centre
165  cellCtrs[own[facei]] += pyr3Vol*pc;
166 
167  // Accumulate face-pyramid volume
168  cellVols[own[facei]] += pyr3Vol;
169  }
170 
171  forAll(nei, facei)
172  {
173  const solveVector fc(fCtrs[facei]);
174  const solveVector fA(fAreas[facei]);
175 
176  // Calculate 3*face-pyramid volume
177  solveScalar pyr3Vol =
178  fA & (cEst[nei[facei]] - fc);
179 
180  // Calculate face-pyramid centre
181  solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[nei[facei]];
182 
183  // Accumulate volume-weighted face-pyramid centre
184  cellCtrs[nei[facei]] += pyr3Vol*pc;
185 
186  // Accumulate face-pyramid volume
187  cellVols[nei[facei]] += pyr3Vol;
188  }
189 
190  forAll(cellCtrs, celli)
191  {
192  if (mag(cellVols[celli]) > VSMALL)
193  {
194  cellCtrs[celli] /= cellVols[celli];
195  }
196  else
197  {
198  cellCtrs[celli] = cEst[celli];
199  }
200  }
201 
202  cellVols *= (1.0/3.0);
203 }
204 
205 
207 (
208  const primitiveMesh& mesh,
209  const pointField& p,
210  const vectorField& fCtrs,
211  const vectorField& fAreas,
212 
213  const label facei,
214  const point& ownCc,
215  const point& neiCc
216 )
217 {
218  vector Cpf = fCtrs[facei] - ownCc;
219  vector d = neiCc - ownCc;
220 
221  // Skewness vector
222  vector sv =
223  Cpf
224  - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
225  vector svHat = sv/(mag(sv) + ROOTVSMALL);
226 
227  // Normalisation distance calculated as the approximate distance
228  // from the face centre to the edge of the face in the direction
229  // of the skewness
230  scalar fd = 0.2*mag(d) + ROOTVSMALL;
231  const face& f = mesh.faces()[facei];
232  forAll(f, pi)
233  {
234  fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei])));
235  }
236 
237  // Normalised skewness
238  return mag(sv)/fd;
239 }
240 
241 
243 (
244  const primitiveMesh& mesh,
245  const pointField& p,
246  const vectorField& fCtrs,
247  const vectorField& fAreas,
248 
249  const label facei,
250  const point& ownCc
251 )
252 {
253  vector Cpf = fCtrs[facei] - ownCc;
254 
255  vector normal = normalised(fAreas[facei]);
256  vector d = normal*(normal & Cpf);
257 
258 
259  // Skewness vector
260  vector sv =
261  Cpf
262  - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d;
263  vector svHat = sv/(mag(sv) + ROOTVSMALL);
264 
265  // Normalisation distance calculated as the approximate distance
266  // from the face centre to the edge of the face in the direction
267  // of the skewness
268  scalar fd = 0.4*mag(d) + ROOTVSMALL;
269  const face& f = mesh.faces()[facei];
270  forAll(f, pi)
271  {
272  fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei])));
273  }
274 
275  // Normalised skewness
276  return mag(sv)/fd;
277 }
278 
279 
281 (
282  const point& ownCc,
283  const point& neiCc,
284  const vector& s
285 )
286 {
287  vector d = neiCc - ownCc;
288 
289  return (d & s)/(mag(d)*mag(s) + ROOTVSMALL);
290 }
291 
292 
293 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
294 
296 (
297  const primitiveMesh& mesh,
298  const vectorField& areas,
299  const vectorField& cc
300 )
301 {
302  const labelList& own = mesh.faceOwner();
303  const labelList& nei = mesh.faceNeighbour();
304 
305  tmp<scalarField> tortho(new scalarField(mesh.nInternalFaces()));
306  scalarField& ortho = tortho.ref();
307 
308  // Internal faces
309  forAll(nei, facei)
310  {
311  ortho[facei] = faceOrthogonality
312  (
313  cc[own[facei]],
314  cc[nei[facei]],
315  areas[facei]
316  );
317  }
318 
319  return tortho;
320 }
321 
322 
324 (
325  const primitiveMesh& mesh,
326  const pointField& p,
327  const vectorField& fCtrs,
328  const vectorField& fAreas,
329  const vectorField& cellCtrs
330 )
331 {
332  const labelList& own = mesh.faceOwner();
333  const labelList& nei = mesh.faceNeighbour();
334 
335  tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
336  scalarField& skew = tskew.ref();
337 
338  forAll(nei, facei)
339  {
340  skew[facei] = faceSkewness
341  (
342  mesh,
343  p,
344  fCtrs,
345  fAreas,
346 
347  facei,
348  cellCtrs[own[facei]],
349  cellCtrs[nei[facei]]
350  );
351  }
352 
353 
354  // Boundary faces: consider them to have only skewness error.
355  // (i.e. treat as if mirror cell on other side)
356 
357  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
358  {
359  skew[facei] = boundaryFaceSkewness
360  (
361  mesh,
362  p,
363  fCtrs,
364  fAreas,
365  facei,
366  cellCtrs[own[facei]]
367  );
368  }
369 
370  return tskew;
371 }
372 
373 
375 (
376  const primitiveMesh& mesh,
377  const pointField& points,
378  const vectorField& ctrs,
379 
380  scalarField& ownPyrVol,
381  scalarField& neiPyrVol
382 )
383 {
384  const labelList& own = mesh.faceOwner();
385  const labelList& nei = mesh.faceNeighbour();
386  const faceList& f = mesh.faces();
387 
388  ownPyrVol.setSize(mesh.nFaces());
389  neiPyrVol.setSize(mesh.nInternalFaces());
390 
391  forAll(f, facei)
392  {
393  // Create the owner pyramid
394  ownPyrVol[facei] = -pyramidPointFaceRef
395  (
396  f[facei],
397  ctrs[own[facei]]
398  ).mag(points);
399 
400  if (mesh.isInternalFace(facei))
401  {
402  // Create the neighbour pyramid - it will have positive volume
403  neiPyrVol[facei] = pyramidPointFaceRef
404  (
405  f[facei],
406  ctrs[nei[facei]]
407  ).mag(points);
408  }
409  }
410 }
411 
412 
414 (
415  const primitiveMesh& mesh,
416  const Vector<label>& meshD,
417  const vectorField& areas,
418  const scalarField& vols,
419 
420  scalarField& openness,
421  scalarField& aratio
422 )
423 {
424  const labelList& own = mesh.faceOwner();
425  const labelList& nei = mesh.faceNeighbour();
426 
427  // Loop through cell faces and sum up the face area vectors for each cell.
428  // This should be zero in all vector components
429 
430  vectorField sumClosed(mesh.nCells(), Zero);
431  vectorField sumMagClosed(mesh.nCells(), Zero);
432 
433  forAll(own, facei)
434  {
435  // Add to owner
436  sumClosed[own[facei]] += areas[facei];
437  sumMagClosed[own[facei]] += cmptMag(areas[facei]);
438  }
439 
440  forAll(nei, facei)
441  {
442  // Subtract from neighbour
443  sumClosed[nei[facei]] -= areas[facei];
444  sumMagClosed[nei[facei]] += cmptMag(areas[facei]);
445  }
446 
447 
448  label nDims = 0;
449  for (direction dir = 0; dir < vector::nComponents; dir++)
450  {
451  if (meshD[dir] == 1)
452  {
453  nDims++;
454  }
455  }
456 
457 
458  // Check the sums
459  openness.setSize(mesh.nCells());
460  aratio.setSize(mesh.nCells());
461 
462  forAll(sumClosed, celli)
463  {
464  scalar maxOpenness = 0;
465 
466  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
467  {
468  maxOpenness = max
469  (
470  maxOpenness,
471  mag(sumClosed[celli][cmpt])
472  /(sumMagClosed[celli][cmpt] + ROOTVSMALL)
473  );
474  }
475  openness[celli] = maxOpenness;
476 
477  // Calculate the aspect ration as the maximum of Cartesian component
478  // aspect ratio to the total area hydraulic area aspect ratio
479  scalar minCmpt = VGREAT;
480  scalar maxCmpt = -VGREAT;
481  for (direction dir = 0; dir < vector::nComponents; dir++)
482  {
483  if (meshD[dir] == 1)
484  {
485  minCmpt = min(minCmpt, sumMagClosed[celli][dir]);
486  maxCmpt = max(maxCmpt, sumMagClosed[celli][dir]);
487  }
488  }
489 
490  scalar aspectRatio = maxCmpt/(minCmpt + ROOTVSMALL);
491  if (nDims == 3)
492  {
493  scalar v = max(ROOTVSMALL, vols[celli]);
494 
495  aspectRatio = max
496  (
497  aspectRatio,
498  1.0/6.0*cmptSum(sumMagClosed[celli])/pow(v, 2.0/3.0)
499  );
500  }
501 
502  aratio[celli] = aspectRatio;
503  }
504 }
505 
506 
508 (
509  const scalar maxSin,
510  const primitiveMesh& mesh,
511  const pointField& p,
512  const vectorField& faceAreas
513 )
514 {
515  const faceList& fcs = mesh.faces();
516 
517  vectorField faceNormals(faceAreas);
518  faceNormals /= mag(faceNormals) + ROOTVSMALL;
519 
520  tmp<scalarField> tfaceAngles(new scalarField(mesh.nFaces()));
521  scalarField& faceAngles = tfaceAngles.ref();
522 
523 
524  forAll(fcs, facei)
525  {
526  const face& f = fcs[facei];
527 
528  // Get edge from f[0] to f[size-1];
529  vector ePrev(p[f.first()] - p[f.last()]);
530  scalar magEPrev = mag(ePrev);
531  ePrev /= magEPrev + ROOTVSMALL;
532 
533  scalar maxEdgeSin = 0.0;
534 
535  forAll(f, fp0)
536  {
537  // Get vertex after fp
538  label fp1 = f.fcIndex(fp0);
539 
540  // Normalized vector between two consecutive points
541  vector e10(p[f[fp1]] - p[f[fp0]]);
542  scalar magE10 = mag(e10);
543  e10 /= magE10 + ROOTVSMALL;
544 
545  if (magEPrev > SMALL && magE10 > SMALL)
546  {
547  vector edgeNormal = ePrev ^ e10;
548  scalar magEdgeNormal = mag(edgeNormal);
549 
550  if (magEdgeNormal < maxSin)
551  {
552  // Edges (almost) aligned -> face is ok.
553  }
554  else
555  {
556  // Check normal
557  edgeNormal /= magEdgeNormal;
558 
559  if ((edgeNormal & faceNormals[facei]) < SMALL)
560  {
561  maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
562  }
563  }
564  }
565 
566  ePrev = e10;
567  magEPrev = magE10;
568  }
569 
570  faceAngles[facei] = maxEdgeSin;
571  }
572 
573  return tfaceAngles;
574 }
575 
576 
578 (
579  const primitiveMesh& mesh,
580  const pointField& p,
581  const vectorField& fCtrs,
582  const vectorField& faceAreas
583 )
584 {
585  const faceList& fcs = mesh.faces();
586 
587  // Areas are calculated as the sum of areas. (see
588  // primitiveMeshFaceCentresAndAreas.C)
589  scalarField magAreas(mag(faceAreas));
590 
591  tmp<scalarField> tfaceFlatness(new scalarField(mesh.nFaces(), 1.0));
592  scalarField& faceFlatness = tfaceFlatness.ref();
593 
594  typedef Vector<solveScalar> solveVector;
595 
596  forAll(fcs, facei)
597  {
598  const face& f = fcs[facei];
599 
600  if (f.size() > 3 && magAreas[facei] > ROOTVSMALL)
601  {
602  const solveVector fc = fCtrs[facei];
603 
604  // Calculate the sum of magnitude of areas and compare to magnitude
605  // of sum of areas.
606 
607  solveScalar sumA = 0.0;
608 
609  forAll(f, fp)
610  {
611  const solveVector thisPoint = p[f[fp]];
612  const solveVector nextPoint = p[f.nextLabel(fp)];
613 
614  // Triangle around fc.
615  solveVector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
616  sumA += mag(n);
617  }
618 
619  faceFlatness[facei] = magAreas[facei]/(sumA + ROOTVSMALL);
620  }
621  }
622 
623  return tfaceFlatness;
624 }
625 
626 
628 (
629  const primitiveMesh& mesh,
630  const Vector<label>& meshD,
631  const vectorField& faceAreas,
632  const bitSet& internalOrCoupledFace
633 )
634 {
635  // Determine number of dimensions and (for 2D) missing dimension
636  label nDims = 0;
637  label twoD = -1;
638  for (direction dir = 0; dir < vector::nComponents; dir++)
639  {
640  if (meshD[dir] == 1)
641  {
642  nDims++;
643  }
644  else
645  {
646  twoD = dir;
647  }
648  }
649 
650  tmp<scalarField> tcellDeterminant(new scalarField(mesh.nCells()));
651  scalarField& cellDeterminant = tcellDeterminant.ref();
652 
653  const cellList& c = mesh.cells();
654 
655  if (nDims == 1)
656  {
657  cellDeterminant = 1.0;
658  }
659  else
660  {
661  forAll(c, celli)
662  {
663  const labelList& curFaces = c[celli];
664 
665  // Calculate local normalization factor
666  scalar avgArea = 0;
667 
668  label nInternalFaces = 0;
669 
670  forAll(curFaces, i)
671  {
672  if (internalOrCoupledFace.test(curFaces[i]))
673  {
674  avgArea += mag(faceAreas[curFaces[i]]);
675 
676  nInternalFaces++;
677  }
678  }
679 
680  if (nInternalFaces == 0 || avgArea < ROOTVSMALL)
681  {
682  cellDeterminant[celli] = 0;
683  }
684  else
685  {
686  avgArea /= nInternalFaces;
687 
688  symmTensor areaTensor(Zero);
689 
690  forAll(curFaces, i)
691  {
692  if (internalOrCoupledFace.test(curFaces[i]))
693  {
694  areaTensor += sqr(faceAreas[curFaces[i]]/avgArea);
695  }
696  }
697 
698  if (nDims == 2)
699  {
700  // Add the missing eigenvector (such that it does not
701  // affect the determinant)
702  if (twoD == 0)
703  {
704  areaTensor.xx() = 1;
705  }
706  else if (twoD == 1)
707  {
708  areaTensor.yy() = 1;
709  }
710  else
711  {
712  areaTensor.zz() = 1;
713  }
714  }
715 
716  // Note:
717  // - normalise to be 0..1 (since cube has eigenvalues 2 2 2)
718  // - we use the determinant (i.e. 3rd invariant) and not e.g.
719  // condition number (= max ev / min ev) since we are
720  // interested in the minimum connectivity and not the
721  // uniformity. Using the condition number on corner cells
722  // leads to uniformity 1 i.e. equally bad in all three
723  // directions which is not what we want.
724  cellDeterminant[celli] = mag(det(areaTensor))/8.0;
725  }
726  }
727  }
728 
729  return tcellDeterminant;
730 }
731 
732 
733 // ************************************************************************* //
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::primitiveMeshTools::cellDeterminant
static tmp< scalarField > cellDeterminant(const primitiveMesh &mesh, const Vector< label > &directions, const vectorField &faceAreas, const bitSet &internalOrCoupledFace)
Generate cell determinant field. Normalised to 1 for an internal cube.
Definition: primitiveMeshTools.C:628
Foam::SymmTensor< scalar >
p
volScalarField & p
Definition: createFieldRefs.H:8
s
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))
Definition: gmvOutputSpray.H:25
Foam::bitSet
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:63
Foam::skew
dimensionedTensor skew(const dimensionedTensor &dt)
Definition: dimensionedTensor.C:138
primitiveMeshTools.H
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::primitiveMeshTools::faceConcavity
static tmp< scalarField > faceConcavity(const scalar maxSin, const primitiveMesh &mesh, const pointField &p, const vectorField &faceAreas)
Generate face concavity field. Returns per face the (sin of the)
Definition: primitiveMeshTools.C:508
PrecisionAdaptor.H
Foam::primitiveMeshTools::faceOrthogonality
static tmp< scalarField > faceOrthogonality(const primitiveMesh &mesh, const vectorField &fAreas, const vectorField &cellCtrs)
Generate non-orthogonality field (internal faces only)
Definition: primitiveMeshTools.C:296
primitiveMesh.H
Foam::bitSet::test
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
Definition: bitSetI.H:520
syncTools.H
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
Foam::SymmTensor::xx
const Cmpt & xx() const
Definition: SymmTensorI.H:84
Foam::cmptMag
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:400
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
Foam::primitiveMeshTools::makeCellCentresAndVols
static void makeCellCentresAndVols(const primitiveMesh &mesh, const vectorField &fCtrs, const vectorField &fAreas, vectorField &cellCtrs, scalarField &cellVols)
Calculate cell centres and volumes from face properties.
Definition: primitiveMeshTools.C:107
faceNormals
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf())
Foam::Field< vector >
cellVols
const scalarField & cellVols
Definition: temperatureAndPressureVariables.H:51
Foam::primitiveMeshTools::faceFlatness
static tmp< scalarField > faceFlatness(const primitiveMesh &mesh, const pointField &p, const vectorField &fCtrs, const vectorField &faceAreas)
Generate face flatness field. Compares the individual triangles'.
Definition: primitiveMeshTools.C:578
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:75
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::cmptSum
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
Return the sum of components of a SphericalTensor.
Definition: SphericalTensorI.H:164
Foam::normalised
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:487
Foam::primitiveMeshTools::faceSkewness
static tmp< scalarField > faceSkewness(const primitiveMesh &mesh, const pointField &points, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs)
Generate skewness field.
Definition: primitiveMeshTools.C:324
Foam::PrecisionAdaptor
A non-const Field/List wrapper with possible data conversion.
Definition: PrecisionAdaptor.H:186
Foam::primitiveMeshTools::makeFaceCentresAndAreas
static void makeFaceCentresAndAreas(const primitiveMesh &mesh, const pointField &p, vectorField &fCtrs, vectorField &fAreas)
Calculate face centres and areas.
Definition: primitiveMeshTools.C:38
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::SymmTensor::yy
const Cmpt & yy() const
Definition: SymmTensorI.H:108
f
labelList f(nPoints)
Foam::Vector
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
Definition: Vector.H:62
Foam::List< face >
Foam::primitiveMeshTools::cellClosedness
static void cellClosedness(const primitiveMesh &mesh, const Vector< label > &meshD, const vectorField &areas, const scalarField &vols, scalarField &openness, scalarField &aratio)
Generate cell openness and cell ascpect ratio field.
Definition: primitiveMeshTools.C:414
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::primitiveMeshTools::boundaryFaceSkewness
static scalar boundaryFaceSkewness(const primitiveMesh &mesh, const pointField &p, const vectorField &fCtrs, const vectorField &fAreas, const label facei, const point &ownCc)
Skewness of single boundary face.
Definition: primitiveMeshTools.C:243
Foam::direction
uint8_t direction
Definition: direction.H:52
Foam::SymmTensor::zz
const Cmpt & zz() const
Definition: SymmTensorI.H:132
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:72
Foam::det
dimensionedScalar det(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:62
Foam::refPtr< Container< Type > >::ref
Container< Type > & ref() const
Definition: refPtrI.H:203
Foam::pyramidPointFaceRef
pyramid< point, const point &, const face & > pyramidPointFaceRef
A pyramid using referred points and faces.
Definition: pyramidPointFaceRef.H:48
pyramidPointFaceRef.H
Foam::primitiveMeshTools::facePyramidVolume
static void facePyramidVolume(const primitiveMesh &mesh, const pointField &points, const vectorField &cellCtrs, scalarField &ownPyrVol, scalarField &neiPyrVol)
Generate face pyramid volume fields.
Definition: primitiveMeshTools.C:375
Foam::primitiveMesh
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:78