polyMeshGeometry.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) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2015-2019 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 "polyMeshGeometry.H"
31 #include "pyramidPointFaceRef.H"
32 #include "tetPointRef.H"
33 #include "syncTools.H"
34 #include "unitConversion.H"
35 #include "primitiveMeshTools.H"
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41  defineTypeNameAndDebug(polyMeshGeometry, 0);
42 }
43 
44 
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 
47 void Foam::polyMeshGeometry::updateFaceCentresAndAreas
48 (
49  const pointField& p,
50  const labelList& changedFaces
51 )
52 {
53  const faceList& fs = mesh_.faces();
54 
55  for (const label facei : changedFaces)
56  {
57  const labelList& f = fs[facei];
58  label nPoints = f.size();
59 
60  // If the face is a triangle, do a direct calculation for efficiency
61  // and to avoid round-off error-related problems
62  if (nPoints == 3)
63  {
64  faceCentres_[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
65  faceAreas_[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
66  }
67  else
68  {
69  vector sumN = Zero;
70  scalar sumA = Zero;
71  vector sumAc = Zero;
72 
73  point fCentre = p[f[0]];
74  for (label pi = 1; pi < nPoints; ++pi)
75  {
76  fCentre += p[f[pi]];
77  }
78 
79  fCentre /= nPoints;
80 
81  for (label pi = 0; pi < nPoints; ++pi)
82  {
83  const point& nextPoint = p[f[(pi + 1) % nPoints]];
84 
85  vector c = p[f[pi]] + nextPoint + fCentre;
86  vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
87  scalar a = mag(n);
88 
89  sumN += n;
90  sumA += a;
91  sumAc += a*c;
92  }
93 
94  faceCentres_[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
95  faceAreas_[facei] = 0.5*sumN;
96  }
97  }
98 }
99 
100 
101 void Foam::polyMeshGeometry::updateCellCentresAndVols
102 (
103  const labelList& changedCells,
104  const labelList& changedFaces
105 )
106 {
107  const labelList& own = mesh().faceOwner();
108  const cellList& cells = mesh().cells();
109 
110  // Clear the fields for accumulation
111  UIndirectList<vector>(cellCentres_, changedCells) = Zero;
112  UIndirectList<scalar>(cellVolumes_, changedCells) = Zero;
113 
114 
115  // Re-calculate the changed cell centres and volumes
116  for (const label celli : changedCells)
117  {
118  const labelList& cFaces = cells[celli];
119 
120  // Estimate the cell centre and bounding box using the face centres
121  vector cEst(Zero);
122  boundBox bb(boundBox::invertedBox);
123 
124  for (const label facei : cFaces)
125  {
126  const point& fc = faceCentres_[facei];
127  cEst += fc;
128  bb.add(fc);
129  }
130  cEst /= cFaces.size();
131 
132 
133  // Sum up the face-pyramid contributions
134  for (const label facei : cFaces)
135  {
136  // Calculate 3* the face-pyramid volume
137  scalar pyr3Vol = faceAreas_[facei] & (faceCentres_[facei] - cEst);
138 
139  if (own[facei] != celli)
140  {
141  pyr3Vol = -pyr3Vol;
142  }
143 
144  // Accumulate face-pyramid volume
145  cellVolumes_[celli] += pyr3Vol;
146 
147  // Calculate the face-pyramid centre
148  const vector pCtr = (3.0/4.0)*faceCentres_[facei] + (1.0/4.0)*cEst;
149 
150  // Accumulate volume-weighted face-pyramid centre
151  cellCentres_[celli] += pyr3Vol*pCtr;
152  }
153 
154  // Average the accumulated quantities
155 
156  if (mag(cellVolumes_[celli]) > VSMALL)
157  {
158  point cc = cellCentres_[celli] / cellVolumes_[celli];
159 
160  // Do additional check for collapsed cells since some volumes
161  // (e.g. 1e-33) do not trigger above but do return completely
162  // wrong cell centre
163  if (bb.contains(cc))
164  {
165  cellCentres_[celli] = cc;
166  }
167  else
168  {
169  cellCentres_[celli] = cEst;
170  }
171  }
172  else
173  {
174  cellCentres_[celli] = cEst;
175  }
176 
177  cellVolumes_[celli] *= (1.0/3.0);
178  }
179 }
180 
181 
183 (
184  const polyMesh& mesh,
185  const labelList& changedFaces
186 )
187 {
188  const labelList& own = mesh.faceOwner();
189  const labelList& nei = mesh.faceNeighbour();
190 
191  labelHashSet affectedCells(2*changedFaces.size());
192 
193  for (const label facei : changedFaces)
194  {
195  affectedCells.insert(own[facei]);
196 
197  if (mesh.isInternalFace(facei))
198  {
199  affectedCells.insert(nei[facei]);
200  }
201  }
202  return affectedCells.toc();
203 }
204 
205 
206 Foam::scalar Foam::polyMeshGeometry::checkNonOrtho
207 (
208  const polyMesh& mesh,
209  const bool report,
210  const scalar severeNonorthogonalityThreshold,
211  const label facei,
212  const vector& s, // face area vector
213  const vector& d, // cc-cc vector
214 
215  label& severeNonOrth,
216  label& errorNonOrth,
217  labelHashSet* setPtr
218 )
219 {
220  scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
221 
222  if (dDotS < severeNonorthogonalityThreshold)
223  {
224  label nei = -1;
225 
226  if (mesh.isInternalFace(facei))
227  {
228  nei = mesh.faceNeighbour()[facei];
229  }
230 
231  if (dDotS > SMALL)
232  {
233  if (report)
234  {
235  // Severe non-orthogonality but mesh still OK
236  Pout<< "Severe non-orthogonality for face " << facei
237  << " between cells " << mesh.faceOwner()[facei]
238  << " and " << nei
239  << ": Angle = "
240  << radToDeg(::acos(dDotS))
241  << " deg." << endl;
242  }
243 
244  ++severeNonOrth;
245  }
246  else
247  {
248  // Non-orthogonality greater than 90 deg
249  if (report)
250  {
252  << "Severe non-orthogonality detected for face "
253  << facei
254  << " between cells " << mesh.faceOwner()[facei]
255  << " and " << nei
256  << ": Angle = "
257  << radToDeg(::acos(dDotS))
258  << " deg." << endl;
259  }
260 
261  ++errorNonOrth;
262  }
263 
264  if (setPtr)
265  {
266  setPtr->insert(facei);
267  }
268  }
269  return dDotS;
270 }
271 
272 
273 // Create the neighbour pyramid - it will have positive volume
274 bool Foam::polyMeshGeometry::checkFaceTet
275 (
276  const polyMesh& mesh,
277  const bool report,
278  const scalar minTetQuality,
279  const pointField& p,
280  const label facei,
281  const point& fc, // face centre
282  const point& cc, // cell centre
283 
284  labelHashSet* setPtr
285 )
286 {
287  const face& f = mesh.faces()[facei];
288 
289  forAll(f, fp)
290  {
291  scalar tetQual = tetPointRef
292  (
293  p[f[fp]],
294  p[f.nextLabel(fp)],
295  fc,
296  cc
297  ).quality();
298 
299  if (tetQual < minTetQuality)
300  {
301  if (report)
302  {
303  Pout<< "bool polyMeshGeometry::checkFaceTets("
304  << "const bool, const scalar, const pointField&"
305  << ", const pointField&"
306  << ", const labelList&, labelHashSet*) : "
307  << "face " << facei
308  << " has a triangle that points the wrong way." << nl
309  << "Tet quality: " << tetQual
310  << " Face " << facei
311  << endl;
312  }
313  if (setPtr)
314  {
315  setPtr->insert(facei);
316  }
317  return true;
318  }
319  }
320 
321  return false;
322 }
323 
324 
325 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
326 
328 :
329  mesh_(mesh)
330 {
331  correct();
332 }
333 
334 
335 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
336 
338 {
339  faceAreas_ = mesh_.faceAreas();
340  faceCentres_ = mesh_.faceCentres();
341  cellCentres_ = mesh_.cellCentres();
342  cellVolumes_ = mesh_.cellVolumes();
343 }
344 
345 
347 (
348  const pointField& p,
349  const labelList& changedFaces
350 )
351 {
352  // Update face quantities
353  updateFaceCentresAndAreas(p, changedFaces);
354  // Update cell quantities from face quantities
355  updateCellCentresAndVols(affectedCells(mesh_, changedFaces), changedFaces);
356 }
357 
358 
360 (
361  const bool report,
362  const scalar orthWarn,
363  const polyMesh& mesh,
364  const vectorField& cellCentres,
365  const vectorField& faceAreas,
366  const labelList& checkFaces,
367  const List<labelPair>& baffles,
368  labelHashSet* setPtr
369 )
370 {
371  // for all internal and coupled faces check that the d dot S product
372  // is positive
373 
374  const labelList& own = mesh.faceOwner();
375  const labelList& nei = mesh.faceNeighbour();
377 
378  // Severe nonorthogonality threshold
379  const scalar severeNonorthogonalityThreshold = ::cos(degToRad(orthWarn));
380 
381  // Calculate coupled cell centre
382  pointField neiCc(mesh.nBoundaryFaces());
383 
384  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
385  {
386  neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
387  }
388 
390 
391  scalar minDDotS = GREAT;
392 
393  scalar sumDDotS = 0;
394  label nDDotS = 0;
395 
396  label severeNonOrth = 0;
397 
398  label errorNonOrth = 0;
399 
400  for (const label facei : checkFaces)
401  {
402  const point& ownCc = cellCentres[own[facei]];
403 
404  if (mesh.isInternalFace(facei))
405  {
406  scalar dDotS = checkNonOrtho
407  (
408  mesh,
409  report,
410  severeNonorthogonalityThreshold,
411  facei,
412  faceAreas[facei],
413  cellCentres[nei[facei]] - ownCc,
414 
415  severeNonOrth,
416  errorNonOrth,
417  setPtr
418  );
419 
420  if (dDotS < minDDotS)
421  {
422  minDDotS = dDotS;
423  }
424 
425  sumDDotS += dDotS;
426  ++nDDotS;
427  }
428  else
429  {
430  const label patchi = patches.whichPatch(facei);
431 
432  if (patches[patchi].coupled())
433  {
434  scalar dDotS = checkNonOrtho
435  (
436  mesh,
437  report,
438  severeNonorthogonalityThreshold,
439  facei,
440  faceAreas[facei],
441  neiCc[facei-mesh.nInternalFaces()] - ownCc,
442 
443  severeNonOrth,
444  errorNonOrth,
445  setPtr
446  );
447 
448  if (dDotS < minDDotS)
449  {
450  minDDotS = dDotS;
451  }
452 
453  sumDDotS += dDotS;
454  ++nDDotS;
455  }
456  }
457  }
458 
459  for (const labelPair& baffle : baffles)
460  {
461  const label face0 = baffle.first();
462  const label face1 = baffle.second();
463 
464  const point& ownCc = cellCentres[own[face0]];
465 
466  scalar dDotS = checkNonOrtho
467  (
468  mesh,
469  report,
470  severeNonorthogonalityThreshold,
471  face0,
472  faceAreas[face0],
473  cellCentres[own[face1]] - ownCc,
474 
475  severeNonOrth,
476  errorNonOrth,
477  setPtr
478  );
479 
480  if (dDotS < minDDotS)
481  {
482  minDDotS = dDotS;
483  }
484 
485  sumDDotS += dDotS;
486  ++nDDotS;
487  }
488 
489  reduce(minDDotS, minOp<scalar>());
490  reduce(sumDDotS, sumOp<scalar>());
491  reduce(nDDotS, sumOp<label>());
492  reduce(severeNonOrth, sumOp<label>());
493  reduce(errorNonOrth, sumOp<label>());
494 
495  // Only report if there are some internal faces
496  if (nDDotS > 0)
497  {
498  if (report && minDDotS < severeNonorthogonalityThreshold)
499  {
500  Info<< "Number of non-orthogonality errors: " << errorNonOrth
501  << ". Number of severely non-orthogonal faces: "
502  << severeNonOrth << "." << endl;
503  }
504  }
505 
506  if (report)
507  {
508  if (nDDotS > 0)
509  {
510  Info<< "Mesh non-orthogonality Max: "
511  << radToDeg(::acos(minDDotS))
512  << " average: " << radToDeg(::acos(sumDDotS/nDDotS))
513  << endl;
514  }
515  }
516 
517  if (errorNonOrth > 0)
518  {
519  if (report)
520  {
522  << "Error in non-orthogonality detected" << endl;
523  }
524 
525  return true;
526  }
527 
528  if (report)
529  {
530  Info<< "Non-orthogonality check OK.\n" << endl;
531  }
532 
533  return false;
534 }
535 
536 
538 (
539  const bool report,
540  const scalar minPyrVol,
541  const polyMesh& mesh,
542  const vectorField& cellCentres,
543  const pointField& p,
544  const labelList& checkFaces,
545  const List<labelPair>& baffles,
546  labelHashSet* setPtr
547 )
548 {
549  // check whether face area vector points to the cell with higher label
550  const labelList& own = mesh.faceOwner();
551  const labelList& nei = mesh.faceNeighbour();
552 
553  const faceList& f = mesh.faces();
554 
555  label nErrorPyrs = 0;
556 
557  for (const label facei : checkFaces)
558  {
559  // Create the owner pyramid - it will have negative volume
560  scalar pyrVol = pyramidPointFaceRef
561  (
562  f[facei],
563  cellCentres[own[facei]]
564  ).mag(p);
565 
566  if (pyrVol > -minPyrVol)
567  {
568  ++nErrorPyrs;
569 
570  if (report)
571  {
572  Pout<< "bool polyMeshGeometry::checkFacePyramids("
573  << "const bool, const scalar, const pointField&"
574  << ", const labelList&, labelHashSet*): "
575  << "face " << facei << " points the wrong way. " << endl
576  << "Pyramid volume: " << -pyrVol
577  << " Face " << f[facei] << " area: " << f[facei].mag(p)
578  << " Owner cell: " << own[facei] << endl
579  << "Owner cell vertex labels: "
580  << mesh.cells()[own[facei]].labels(f)
581  << endl;
582  }
583  if (setPtr)
584  {
585  setPtr->insert(facei);
586  }
587  }
588 
589  if (mesh.isInternalFace(facei))
590  {
591  // Create the neighbour pyramid - it will have positive volume
592  scalar pyrVol =
593  pyramidPointFaceRef(f[facei], cellCentres[nei[facei]]).mag(p);
594 
595  if (pyrVol < minPyrVol)
596  {
597  ++nErrorPyrs;
598 
599  if (report)
600  {
601  Pout<< "bool polyMeshGeometry::checkFacePyramids("
602  << "const bool, const scalar, const pointField&"
603  << ", const labelList&, labelHashSet*): "
604  << "face " << facei << " points the wrong way. " << endl
605  << "Pyramid volume: " << -pyrVol
606  << " Face " << f[facei] << " area: " << f[facei].mag(p)
607  << " Neighbour cell: " << nei[facei] << endl
608  << "Neighbour cell vertex labels: "
609  << mesh.cells()[nei[facei]].labels(f)
610  << endl;
611  }
612  if (setPtr)
613  {
614  setPtr->insert(facei);
615  }
616  }
617  }
618  }
619 
620  for (const labelPair& baffle : baffles)
621  {
622  const label face0 = baffle.first();
623  const label face1 = baffle.second();
624 
625  const point& ownCc = cellCentres[own[face0]];
626 
627  // Create the owner pyramid - it will have negative volume
628  scalar pyrVolOwn = pyramidPointFaceRef
629  (
630  f[face0],
631  ownCc
632  ).mag(p);
633 
634  if (pyrVolOwn > -minPyrVol)
635  {
636  ++nErrorPyrs;
637 
638  if (report)
639  {
640  Pout<< "bool polyMeshGeometry::checkFacePyramids("
641  << "const bool, const scalar, const pointField&"
642  << ", const labelList&, labelHashSet*): "
643  << "face " << face0 << " points the wrong way. " << endl
644  << "Pyramid volume: " << -pyrVolOwn
645  << " Face " << f[face0] << " area: " << f[face0].mag(p)
646  << " Owner cell: " << own[face0] << endl
647  << "Owner cell vertex labels: "
648  << mesh.cells()[own[face0]].labels(f)
649  << endl;
650  }
651  if (setPtr)
652  {
653  setPtr->insert(face0);
654  }
655  }
656 
657  // Create the neighbour pyramid - it will have positive volume
658  scalar pyrVolNbr =
659  pyramidPointFaceRef(f[face0], cellCentres[own[face1]]).mag(p);
660 
661  if (pyrVolNbr < minPyrVol)
662  {
663  ++nErrorPyrs;
664 
665  if (report)
666  {
667  Pout<< "bool polyMeshGeometry::checkFacePyramids("
668  << "const bool, const scalar, const pointField&"
669  << ", const labelList&, labelHashSet*): "
670  << "face " << face0 << " points the wrong way. " << endl
671  << "Pyramid volume: " << -pyrVolNbr
672  << " Face " << f[face0] << " area: " << f[face0].mag(p)
673  << " Neighbour cell: " << own[face1] << endl
674  << "Neighbour cell vertex labels: "
675  << mesh.cells()[own[face1]].labels(f)
676  << endl;
677  }
678  if (setPtr)
679  {
680  setPtr->insert(face0);
681  }
682  }
683  }
684 
685  reduce(nErrorPyrs, sumOp<label>());
686 
687  if (nErrorPyrs > 0)
688  {
689  if (report)
690  {
692  << "Error in face pyramids: faces pointing the wrong way."
693  << endl;
694  }
695 
696  return true;
697  }
698 
699  if (report)
700  {
701  Info<< "Face pyramids OK.\n" << endl;
702  }
703 
704  return false;
705 }
706 
707 
709 (
710  const bool report,
711  const scalar minTetQuality,
712  const polyMesh& mesh,
713  const vectorField& cellCentres,
714  const vectorField& faceCentres,
715  const pointField& p,
716  const labelList& checkFaces,
717  const List<labelPair>& baffles,
718  labelHashSet* setPtr
719 )
720 {
721  // check whether decomposing each cell into tets results in
722  // positive volume, non-flat tets
723  const labelList& own = mesh.faceOwner();
724  const labelList& nei = mesh.faceNeighbour();
726 
727  // Calculate coupled cell centre
728  pointField neiCc(mesh.nBoundaryFaces());
729 
730  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
731  {
732  neiCc[facei - mesh.nInternalFaces()] = cellCentres[own[facei]];
733  }
734 
736 
737  label nErrorTets = 0;
738 
739  for (const label facei : checkFaces)
740  {
741  // Create the owner pyramid - note: exchange cell and face centre
742  // to get positive volume.
743  bool tetError = checkFaceTet
744  (
745  mesh,
746  report,
747  minTetQuality,
748  p,
749  facei,
750  cellCentres[own[facei]], // face centre
751  faceCentres[facei], // cell centre
752  setPtr
753  );
754 
755  if (tetError)
756  {
757  ++nErrorTets;
758  }
759 
760  if (mesh.isInternalFace(facei))
761  {
762  // Create the neighbour tets - they will have positive volume
763  bool tetError = checkFaceTet
764  (
765  mesh,
766  report,
767  minTetQuality,
768  p,
769  facei,
770  faceCentres[facei], // face centre
771  cellCentres[nei[facei]], // cell centre
772  setPtr
773  );
774 
775  if (tetError)
776  {
777  ++nErrorTets;
778  }
779 
780  if
781  (
783  (
784  mesh,
785  facei,
786  minTetQuality,
787  report
788  ) == -1
789  )
790  {
791  ++nErrorTets;
792  if (setPtr)
793  {
794  setPtr->insert(facei);
795  }
796  }
797  }
798  else
799  {
800  label patchi = patches.whichPatch(facei);
801 
802  if (patches[patchi].coupled())
803  {
804  if
805  (
807  (
808  mesh,
809  facei,
810  neiCc[facei - mesh.nInternalFaces()],
811  minTetQuality,
812  report
813  ) == -1
814  )
815  {
816  ++nErrorTets;
817  if (setPtr)
818  {
819  setPtr->insert(facei);
820  }
821  }
822  }
823  else
824  {
825  if
826  (
828  (
829  mesh,
830  facei,
831  minTetQuality,
832  report
833  ) == -1
834  )
835  {
836  ++nErrorTets;
837  if (setPtr)
838  {
839  setPtr->insert(facei);
840  }
841  }
842  }
843  }
844  }
845 
846  for (const labelPair& baffle : baffles)
847  {
848  const label face0 = baffle.first();
849  const label face1 = baffle.second();
850 
851  bool tetError = checkFaceTet
852  (
853  mesh,
854  report,
855  minTetQuality,
856  p,
857  face0,
858  cellCentres[own[face0]], // face centre
859  faceCentres[face0], // cell centre
860  setPtr
861  );
862 
863  if (tetError)
864  {
865  ++nErrorTets;
866  }
867 
868  // Create the neighbour tets - they will have positive volume
869  tetError = checkFaceTet
870  (
871  mesh,
872  report,
873  minTetQuality,
874  p,
875  face0,
876  faceCentres[face0], // face centre
877  cellCentres[own[face1]], // cell centre
878  setPtr
879  );
880 
881  if (tetError)
882  {
883  ++nErrorTets;
884  }
885 
886  if
887  (
889  (
890  mesh,
891  face0,
892  cellCentres[own[face1]],
893  minTetQuality,
894  report
895  ) == -1
896  )
897  {
898  ++nErrorTets;
899  if (setPtr)
900  {
901  setPtr->insert(face0);
902  }
903  }
904  }
905 
906  reduce(nErrorTets, sumOp<label>());
907 
908  if (nErrorTets > 0)
909  {
910  if (report)
911  {
913  << "Error in face decomposition: negative tets."
914  << endl;
915  }
916 
917  return true;
918  }
919 
920  if (report)
921  {
922  Info<< "Face tets OK.\n" << endl;
923  }
924 
925  return false;
926 }
927 
928 
930 (
931  const bool report,
932  const scalar internalSkew,
933  const scalar boundarySkew,
934  const polyMesh& mesh,
935  const pointField& points,
936  const vectorField& cellCentres,
937  const vectorField& faceCentres,
938  const vectorField& faceAreas,
939  const labelList& checkFaces,
940  const List<labelPair>& baffles,
941  labelHashSet* setPtr
942 )
943 {
944  // Warn if the skew correction vector is more than skew times
945  // larger than the face area vector
946 
947  const labelList& own = mesh.faceOwner();
948  const labelList& nei = mesh.faceNeighbour();
950 
951  // Calculate coupled cell centre
952  pointField neiCc;
953  syncTools::swapBoundaryCellPositions(mesh, cellCentres, neiCc);
954 
955  scalar maxSkew = 0;
956 
957  label nWarnSkew = 0;
958 
959  for (const label facei : checkFaces)
960  {
961  if (mesh.isInternalFace(facei))
962  {
963  scalar skewness = primitiveMeshTools::faceSkewness
964  (
965  mesh,
966  points,
967  faceCentres,
968  faceAreas,
969 
970  facei,
971  cellCentres[own[facei]],
972  cellCentres[nei[facei]]
973  );
974 
975  // Check if the skewness vector is greater than the PN vector.
976  // This does not cause trouble but is a good indication of a poor
977  // mesh.
978  if (skewness > internalSkew)
979  {
980  ++nWarnSkew;
981 
982  if (report)
983  {
984  Pout<< "Severe skewness for face " << facei
985  << " skewness = " << skewness << endl;
986  }
987  if (setPtr)
988  {
989  setPtr->insert(facei);
990  }
991  }
992 
993  maxSkew = max(maxSkew, skewness);
994  }
995  else if (patches[patches.whichPatch(facei)].coupled())
996  {
997  scalar skewness = primitiveMeshTools::faceSkewness
998  (
999  mesh,
1000  points,
1001  faceCentres,
1002  faceAreas,
1003 
1004  facei,
1005  cellCentres[own[facei]],
1006  neiCc[facei-mesh.nInternalFaces()]
1007  );
1008 
1009  // Check if the skewness vector is greater than the PN vector.
1010  // This does not cause trouble but is a good indication of a poor
1011  // mesh.
1012  if (skewness > internalSkew)
1013  {
1014  ++nWarnSkew;
1015 
1016  if (report)
1017  {
1018  Pout<< "Severe skewness for coupled face " << facei
1019  << " skewness = " << skewness << endl;
1020  }
1021  if (setPtr)
1022  {
1023  setPtr->insert(facei);
1024  }
1025  }
1026 
1027  maxSkew = max(maxSkew, skewness);
1028  }
1029  else
1030  {
1032  (
1033  mesh,
1034  points,
1035  faceCentres,
1036  faceAreas,
1037 
1038  facei,
1039  cellCentres[own[facei]]
1040  );
1041 
1042 
1043  // Check if the skewness vector is greater than the PN vector.
1044  // This does not cause trouble but is a good indication of a poor
1045  // mesh.
1046  if (skewness > boundarySkew)
1047  {
1048  ++nWarnSkew;
1049 
1050  if (report)
1051  {
1052  Pout<< "Severe skewness for boundary face " << facei
1053  << " skewness = " << skewness << endl;
1054  }
1055  if (setPtr)
1056  {
1057  setPtr->insert(facei);
1058  }
1059  }
1060 
1061  maxSkew = max(maxSkew, skewness);
1062  }
1063  }
1064 
1065  for (const labelPair& baffle : baffles)
1066  {
1067  const label face0 = baffle.first();
1068  const label face1 = baffle.second();
1069 
1070  const point& ownCc = cellCentres[own[face0]];
1071  const point& neiCc = cellCentres[own[face1]];
1072 
1073  scalar skewness = primitiveMeshTools::faceSkewness
1074  (
1075  mesh,
1076  points,
1077  faceCentres,
1078  faceAreas,
1079 
1080  face0,
1081  ownCc,
1082  neiCc
1083  );
1084 
1085  // Check if the skewness vector is greater than the PN vector.
1086  // This does not cause trouble but is a good indication of a poor
1087  // mesh.
1088  if (skewness > internalSkew)
1089  {
1090  ++nWarnSkew;
1091 
1092  if (report)
1093  {
1094  Pout<< "Severe skewness for face " << face0
1095  << " skewness = " << skewness << endl;
1096  }
1097  if (setPtr)
1098  {
1099  setPtr->insert(face0);
1100  }
1101  }
1102 
1103  maxSkew = max(maxSkew, skewness);
1104  }
1105 
1106 
1107  reduce(maxSkew, maxOp<scalar>());
1108  reduce(nWarnSkew, sumOp<label>());
1109 
1110  if (nWarnSkew > 0)
1111  {
1112  if (report)
1113  {
1115  << 100*maxSkew
1116  << " percent.\nThis may impair the quality of the result." << nl
1117  << nWarnSkew << " highly skew faces detected."
1118  << endl;
1119  }
1120 
1121  return true;
1122  }
1123 
1124  if (report)
1125  {
1126  Info<< "Max skewness = " << 100*maxSkew
1127  << " percent. Face skewness OK.\n" << endl;
1128  }
1129 
1130  return false;
1131 }
1132 
1133 
1136  const bool report,
1137  const scalar warnWeight,
1138  const polyMesh& mesh,
1139  const vectorField& cellCentres,
1140  const vectorField& faceCentres,
1141  const vectorField& faceAreas,
1142  const labelList& checkFaces,
1143  const List<labelPair>& baffles,
1144  labelHashSet* setPtr
1145 )
1146 {
1147  // Warn if the delta factor (0..1) is too large.
1148 
1149  const labelList& own = mesh.faceOwner();
1150  const labelList& nei = mesh.faceNeighbour();
1152 
1153  // Calculate coupled cell centre
1154  pointField neiCc(mesh.nBoundaryFaces());
1155 
1156  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
1157  {
1158  neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
1159  }
1161 
1162 
1163  scalar minWeight = GREAT;
1164 
1165  label nWarnWeight = 0;
1166 
1167  for (const label facei : checkFaces)
1168  {
1169  const point& fc = faceCentres[facei];
1170  const vector& fa = faceAreas[facei];
1171 
1172  scalar dOwn = mag(fa & (fc-cellCentres[own[facei]]));
1173 
1174  if (mesh.isInternalFace(facei))
1175  {
1176  scalar dNei = mag(fa & (cellCentres[nei[facei]]-fc));
1177  scalar weight = min(dNei,dOwn)/(dNei+dOwn+VSMALL);
1178 
1179  if (weight < warnWeight)
1180  {
1181  ++nWarnWeight;
1182 
1183  if (report)
1184  {
1185  Pout<< "Small weighting factor for face " << facei
1186  << " weight = " << weight << endl;
1187  }
1188  if (setPtr)
1189  {
1190  setPtr->insert(facei);
1191  }
1192  }
1193 
1194  minWeight = min(minWeight, weight);
1195  }
1196  else
1197  {
1198  label patchi = patches.whichPatch(facei);
1199 
1200  if (patches[patchi].coupled())
1201  {
1202  scalar dNei = mag(fa & (neiCc[facei-mesh.nInternalFaces()]-fc));
1203  scalar weight = min(dNei,dOwn)/(dNei+dOwn+VSMALL);
1204 
1205  if (weight < warnWeight)
1206  {
1207  ++nWarnWeight;
1208 
1209  if (report)
1210  {
1211  Pout<< "Small weighting factor for face " << facei
1212  << " weight = " << weight << endl;
1213  }
1214  if (setPtr)
1215  {
1216  setPtr->insert(facei);
1217  }
1218  }
1219 
1220  minWeight = min(minWeight, weight);
1221  }
1222  }
1223  }
1224 
1225  for (const labelPair& baffle : baffles)
1226  {
1227  const label face0 = baffle.first();
1228  const label face1 = baffle.second();
1229 
1230  const point& ownCc = cellCentres[own[face0]];
1231  const point& fc = faceCentres[face0];
1232  const vector& fa = faceAreas[face0];
1233 
1234  scalar dOwn = mag(fa & (fc-ownCc));
1235  scalar dNei = mag(fa & (cellCentres[own[face1]]-fc));
1236  scalar weight = min(dNei,dOwn)/(dNei+dOwn+VSMALL);
1237 
1238  if (weight < warnWeight)
1239  {
1240  ++nWarnWeight;
1241 
1242  if (report)
1243  {
1244  Pout<< "Small weighting factor for face " << face0
1245  << " weight = " << weight << endl;
1246  }
1247  if (setPtr)
1248  {
1249  setPtr->insert(face0);
1250  }
1251  }
1252 
1253  minWeight = min(minWeight, weight);
1254  }
1255 
1256  reduce(minWeight, minOp<scalar>());
1257  reduce(nWarnWeight, sumOp<label>());
1258 
1259  if (minWeight < warnWeight)
1260  {
1261  if (report)
1262  {
1264  << minWeight << '.' << nl
1265  << nWarnWeight << " faces with small weights detected."
1266  << endl;
1267  }
1268 
1269  return true;
1270  }
1271 
1272  if (report)
1273  {
1274  Info<< "Min weight = " << minWeight
1275  << ". Weights OK.\n" << endl;
1276  }
1277 
1278  return false;
1279 }
1280 
1281 
1284  const bool report,
1285  const scalar warnRatio,
1286  const polyMesh& mesh,
1287  const scalarField& cellVolumes,
1288  const labelList& checkFaces,
1289  const List<labelPair>& baffles,
1290  labelHashSet* setPtr
1291 )
1292 {
1293  // Warn if the volume ratio between neighbouring cells is too large
1294 
1295  const labelList& own = mesh.faceOwner();
1296  const labelList& nei = mesh.faceNeighbour();
1298 
1299  // Calculate coupled cell vol
1300  scalarField neiVols(mesh.nBoundaryFaces());
1301 
1302  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
1303  {
1304  neiVols[facei-mesh.nInternalFaces()] = cellVolumes[own[facei]];
1305  }
1307 
1308 
1309  scalar minRatio = GREAT;
1310 
1311  label nWarnRatio = 0;
1312 
1313  for (const label facei : checkFaces)
1314  {
1315  scalar ownVol = mag(cellVolumes[own[facei]]);
1316 
1317  scalar neiVol = -GREAT;
1318 
1319  if (mesh.isInternalFace(facei))
1320  {
1321  neiVol = mag(cellVolumes[nei[facei]]);
1322  }
1323  else
1324  {
1325  label patchi = patches.whichPatch(facei);
1326 
1327  if (patches[patchi].coupled())
1328  {
1329  neiVol = mag(neiVols[facei-mesh.nInternalFaces()]);
1330  }
1331  }
1332 
1333  if (neiVol >= 0)
1334  {
1335  scalar ratio = min(ownVol, neiVol) / (max(ownVol, neiVol) + VSMALL);
1336 
1337  if (ratio < warnRatio)
1338  {
1339  ++nWarnRatio;
1340 
1341  if (report)
1342  {
1343  Pout<< "Small ratio for face " << facei
1344  << " ratio = " << ratio << endl;
1345  }
1346  if (setPtr)
1347  {
1348  setPtr->insert(facei);
1349  }
1350  }
1351 
1352  minRatio = min(minRatio, ratio);
1353  }
1354  }
1355 
1356  for (const labelPair& baffle : baffles)
1357  {
1358  const label face0 = baffle.first();
1359  const label face1 = baffle.second();
1360 
1361  scalar ownVol = mag(cellVolumes[own[face0]]);
1362 
1363  scalar neiVol = mag(cellVolumes[own[face1]]);
1364 
1365  if (neiVol >= 0)
1366  {
1367  scalar ratio = min(ownVol, neiVol) / (max(ownVol, neiVol) + VSMALL);
1368 
1369  if (ratio < warnRatio)
1370  {
1371  ++nWarnRatio;
1372 
1373  if (report)
1374  {
1375  Pout<< "Small ratio for face " << face0
1376  << " ratio = " << ratio << endl;
1377  }
1378  if (setPtr)
1379  {
1380  setPtr->insert(face0);
1381  }
1382  }
1383 
1384  minRatio = min(minRatio, ratio);
1385  }
1386  }
1387 
1388  reduce(minRatio, minOp<scalar>());
1389  reduce(nWarnRatio, sumOp<label>());
1390 
1391  if (minRatio < warnRatio)
1392  {
1393  if (report)
1394  {
1396  << minRatio << '.' << nl
1397  << nWarnRatio << " faces with small ratios detected."
1398  << endl;
1399  }
1400 
1401  return true;
1402  }
1403 
1404  if (report)
1405  {
1406  Info<< "Min ratio = " << minRatio
1407  << ". Ratios OK.\n" << endl;
1408  }
1409 
1410  return false;
1411 }
1412 
1413 
1414 // Check convexity of angles in a face. Allow a slight non-convexity.
1415 // E.g. maxDeg = 10 allows for angles < 190 (or 10 degrees concavity)
1416 // (if truly concave and points not visible from face centre the face-pyramid
1417 // check in checkMesh will fail)
1420  const bool report,
1421  const scalar maxDeg,
1422  const polyMesh& mesh,
1423  const vectorField& faceAreas,
1424  const pointField& p,
1425  const labelList& checkFaces,
1426  labelHashSet* setPtr
1427 )
1428 {
1429  if (maxDeg < -SMALL || maxDeg > 180+SMALL)
1430  {
1432  << "maxDeg should be [0..180] but is now " << maxDeg
1433  << abort(FatalError);
1434  }
1435 
1436  const scalar maxSin = Foam::sin(degToRad(maxDeg));
1437 
1438  const faceList& fcs = mesh.faces();
1439 
1440  scalar maxEdgeSin = 0.0;
1441 
1442  label nConcave = 0;
1443 
1444  label errorFacei = -1;
1445 
1446  for (const label facei : checkFaces)
1447  {
1448  const face& f = fcs[facei];
1449 
1450  const vector faceNormal = normalised(faceAreas[facei]);
1451 
1452  // Get edge from f[0] to f[size-1];
1453  vector ePrev(p[f.first()] - p[f.last()]);
1454  scalar magEPrev = mag(ePrev);
1455  ePrev /= magEPrev + VSMALL;
1456 
1457  forAll(f, fp0)
1458  {
1459  // Get vertex after fp
1460  label fp1 = f.fcIndex(fp0);
1461 
1462  // Normalized vector between two consecutive points
1463  vector e10(p[f[fp1]] - p[f[fp0]]);
1464  scalar magE10 = mag(e10);
1465  e10 /= magE10 + VSMALL;
1466 
1467  if (magEPrev > SMALL && magE10 > SMALL)
1468  {
1469  vector edgeNormal = ePrev ^ e10;
1470  scalar magEdgeNormal = mag(edgeNormal);
1471 
1472  if (magEdgeNormal < maxSin)
1473  {
1474  // Edges (almost) aligned -> face is ok.
1475  }
1476  else
1477  {
1478  // Check normal
1479  edgeNormal /= magEdgeNormal;
1480 
1481  if ((edgeNormal & faceNormal) < SMALL)
1482  {
1483  if (facei != errorFacei)
1484  {
1485  // Count only one error per face.
1486  errorFacei = facei;
1487  ++nConcave;
1488  }
1489 
1490  maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
1491 
1492  if (setPtr)
1493  {
1494  setPtr->insert(facei);
1495  }
1496  }
1497  }
1498  }
1499 
1500  ePrev = e10;
1501  magEPrev = magE10;
1502  }
1503  }
1504 
1505  reduce(nConcave, sumOp<label>());
1506  reduce(maxEdgeSin, maxOp<scalar>());
1507 
1508  if (report)
1509  {
1510  if (maxEdgeSin > SMALL)
1511  {
1512  scalar maxConcaveDegr =
1513  radToDeg(Foam::asin(Foam::min(1.0, maxEdgeSin)));
1514 
1515  Info<< "There are " << nConcave
1516  << " faces with concave angles between consecutive"
1517  << " edges. Max concave angle = " << maxConcaveDegr
1518  << " degrees.\n" << endl;
1519  }
1520  else
1521  {
1522  Info<< "All angles in faces are convex or less than " << maxDeg
1523  << " degrees concave.\n" << endl;
1524  }
1525  }
1526 
1527  if (nConcave > 0)
1528  {
1529  if (report)
1530  {
1532  << nConcave << " face points with severe concave angle (> "
1533  << maxDeg << " deg) found.\n"
1534  << endl;
1535  }
1536 
1537  return true;
1538  }
1539 
1540  return false;
1541 }
1542 
1543 
1544 // Check twist of faces. Is calculated as the difference between normals of
1545 // individual triangles and the cell-cell centre edge.
1548  const bool report,
1549  const scalar minTwist,
1550  const polyMesh& mesh,
1551  const vectorField& cellCentres,
1552  const vectorField& faceAreas,
1553  const vectorField& faceCentres,
1554  const pointField& p,
1555  const labelList& checkFaces,
1556  labelHashSet* setPtr
1557 )
1558 {
1559  if (minTwist < -1-SMALL || minTwist > 1+SMALL)
1560  {
1562  << "minTwist should be [-1..1] but is now " << minTwist
1563  << abort(FatalError);
1564  }
1565 
1566 
1567  const faceList& fcs = mesh.faces();
1568 
1569  label nWarped = 0;
1570 
1571  const labelList& own = mesh.faceOwner();
1572  const labelList& nei = mesh.faceNeighbour();
1574 
1575  // Calculate coupled cell centre
1576  pointField neiCc(mesh.nBoundaryFaces());
1577 
1578  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
1579  {
1580  neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
1581  }
1583 
1584  for (const label facei : checkFaces)
1585  {
1586  const face& f = fcs[facei];
1587 
1588  if (f.size() > 3)
1589  {
1590  vector nf(Zero);
1591 
1592  if (mesh.isInternalFace(facei))
1593  {
1594  nf =
1595  normalised
1596  (
1597  cellCentres[nei[facei]]
1598  - cellCentres[own[facei]]
1599  );
1600  }
1601  else if (patches[patches.whichPatch(facei)].coupled())
1602  {
1603  nf =
1604  normalised
1605  (
1606  neiCc[facei-mesh.nInternalFaces()]
1607  - cellCentres[own[facei]]
1608  );
1609  }
1610  else
1611  {
1612  nf =
1613  normalised
1614  (
1615  faceCentres[facei]
1616  - cellCentres[own[facei]]
1617  );
1618  }
1619 
1620  if (nf != vector::zero)
1621  {
1622  const point& fc = faceCentres[facei];
1623 
1624  forAll(f, fpI)
1625  {
1626  vector triArea
1627  (
1628  triPointRef
1629  (
1630  p[f[fpI]],
1631  p[f.nextLabel(fpI)],
1632  fc
1633  ).areaNormal()
1634  );
1635 
1636  scalar magTri = mag(triArea);
1637 
1638  if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
1639  {
1640  ++nWarped;
1641 
1642  if (setPtr)
1643  {
1644  setPtr->insert(facei);
1645  }
1646 
1647  break;
1648  }
1649  }
1650  }
1651  }
1652  }
1653 
1654  reduce(nWarped, sumOp<label>());
1655 
1656  if (report)
1657  {
1658  if (nWarped> 0)
1659  {
1660  Info<< "There are " << nWarped
1661  << " faces with cosine of the angle"
1662  << " between triangle normal and face normal less than "
1663  << minTwist << nl << endl;
1664  }
1665  else
1666  {
1667  Info<< "All faces are flat in that the cosine of the angle"
1668  << " between triangle normal and face normal less than "
1669  << minTwist << nl << endl;
1670  }
1671  }
1672 
1673  if (nWarped > 0)
1674  {
1675  if (report)
1676  {
1678  << nWarped << " faces with severe warpage "
1679  << "(cosine of the angle between triangle normal and "
1680  << "face normal < " << minTwist << ") found.\n"
1681  << endl;
1682  }
1683 
1684  return true;
1685  }
1686 
1687  return false;
1688 }
1689 
1690 
1691 // Like checkFaceTwist but compares normals of consecutive triangles.
1694  const bool report,
1695  const scalar minTwist,
1696  const polyMesh& mesh,
1697  const vectorField& faceAreas,
1698  const vectorField& faceCentres,
1699  const pointField& p,
1700  const labelList& checkFaces,
1701  labelHashSet* setPtr
1702 )
1703 {
1704  if (minTwist < -1-SMALL || minTwist > 1+SMALL)
1705  {
1707  << "minTwist should be [-1..1] but is now " << minTwist
1708  << abort(FatalError);
1709  }
1710 
1711  const faceList& fcs = mesh.faces();
1712 
1713  label nWarped = 0;
1714 
1715  for (const label facei : checkFaces)
1716  {
1717  const face& f = fcs[facei];
1718 
1719  if (f.size() > 3)
1720  {
1721  const point& fc = faceCentres[facei];
1722 
1723  // Find starting triangle (at startFp) with non-zero area
1724  label startFp = -1;
1725  vector prevN;
1726 
1727  forAll(f, fp)
1728  {
1729  prevN = triPointRef
1730  (
1731  p[f[fp]],
1732  p[f.nextLabel(fp)],
1733  fc
1734  ).areaNormal();
1735 
1736  scalar magTri = mag(prevN);
1737 
1738  if (magTri > VSMALL)
1739  {
1740  startFp = fp;
1741  prevN /= magTri;
1742  break;
1743  }
1744  }
1745 
1746  if (startFp != -1)
1747  {
1748  label fp = startFp;
1749 
1750  do
1751  {
1752  fp = f.fcIndex(fp);
1753 
1754  vector triN
1755  (
1756  triPointRef
1757  (
1758  p[f[fp]],
1759  p[f.nextLabel(fp)],
1760  fc
1761  ).areaNormal()
1762  );
1763  scalar magTri = mag(triN);
1764 
1765  if (magTri > VSMALL)
1766  {
1767  triN /= magTri;
1768 
1769  if ((prevN & triN) < minTwist)
1770  {
1771  ++nWarped;
1772 
1773  if (setPtr)
1774  {
1775  setPtr->insert(facei);
1776  }
1777 
1778  break;
1779  }
1780 
1781  prevN = triN;
1782  }
1783  else if (minTwist > 0)
1784  {
1785  ++nWarped;
1786 
1787  if (setPtr)
1788  {
1789  setPtr->insert(facei);
1790  }
1791 
1792  break;
1793  }
1794  }
1795  while (fp != startFp);
1796  }
1797  }
1798  }
1799 
1800 
1801  reduce(nWarped, sumOp<label>());
1802 
1803  if (report)
1804  {
1805  if (nWarped> 0)
1806  {
1807  Info<< "There are " << nWarped
1808  << " faces with cosine of the angle"
1809  << " between consecutive triangle normals less than "
1810  << minTwist << nl << endl;
1811  }
1812  else
1813  {
1814  Info<< "All faces are flat in that the cosine of the angle"
1815  << " between consecutive triangle normals is less than "
1816  << minTwist << nl << endl;
1817  }
1818  }
1819 
1820  if (nWarped > 0)
1821  {
1822  if (report)
1823  {
1825  << nWarped << " faces with severe warpage "
1826  << "(cosine of the angle between consecutive triangle normals"
1827  << " < " << minTwist << ") found.\n"
1828  << endl;
1829  }
1830 
1831  return true;
1832  }
1833 
1834  return false;
1835 }
1836 
1837 
1840  const bool report,
1841  const scalar minFlatness,
1842  const polyMesh& mesh,
1843  const vectorField& faceAreas,
1844  const vectorField& faceCentres,
1845  const pointField& p,
1846  const labelList& checkFaces,
1847  labelHashSet* setPtr
1848 )
1849 {
1850  if (minFlatness < -SMALL || minFlatness > 1+SMALL)
1851  {
1853  << "minFlatness should be [0..1] but is now " << minFlatness
1854  << abort(FatalError);
1855  }
1856 
1857  const faceList& fcs = mesh.faces();
1858 
1859  label nWarped = 0;
1860 
1861  for (const label facei : checkFaces)
1862  {
1863  const face& f = fcs[facei];
1864 
1865  if (f.size() > 3)
1866  {
1867  const point& fc = faceCentres[facei];
1868 
1869  // Sum triangle areas
1870  scalar sumArea = 0.0;
1871 
1872  forAll(f, fp)
1873  {
1874  sumArea += triPointRef
1875  (
1876  p[f[fp]],
1877  p[f.nextLabel(fp)],
1878  fc
1879  ).mag();
1880  }
1881 
1882  if (sumArea/mag(faceAreas[facei]) < minFlatness)
1883  {
1884  ++nWarped;
1885  if (setPtr)
1886  {
1887  setPtr->insert(facei);
1888  }
1889  }
1890  }
1891  }
1892 
1893  reduce(nWarped, sumOp<label>());
1894 
1895  if (report)
1896  {
1897  if (nWarped> 0)
1898  {
1899  Info<< "There are " << nWarped
1900  << " faces with area of individual triangles"
1901  << " compared to overall area less than "
1902  << minFlatness << nl << endl;
1903  }
1904  else
1905  {
1906  Info<< "All faces are flat in that the area of individual triangles"
1907  << " compared to overall area is less than "
1908  << minFlatness << nl << endl;
1909  }
1910  }
1911 
1912  if (nWarped > 0)
1913  {
1914  if (report)
1915  {
1917  << nWarped << " non-flat faces "
1918  << "(area of individual triangles"
1919  << " compared to overall area"
1920  << " < " << minFlatness << ") found.\n"
1921  << endl;
1922  }
1923 
1924  return true;
1925  }
1926 
1927  return false;
1928 }
1929 
1930 
1933  const bool report,
1934  const scalar minArea,
1935  const polyMesh& mesh,
1936  const vectorField& faceAreas,
1937  const labelList& checkFaces,
1938  labelHashSet* setPtr
1939 )
1940 {
1941  label nZeroArea = 0;
1942 
1943  for (const label facei : checkFaces)
1944  {
1945  if (mag(faceAreas[facei]) < minArea)
1946  {
1947  ++nZeroArea;
1948  if (setPtr)
1949  {
1950  setPtr->insert(facei);
1951  }
1952  }
1953  }
1954 
1955 
1956  reduce(nZeroArea, sumOp<label>());
1957 
1958  if (report)
1959  {
1960  if (nZeroArea > 0)
1961  {
1962  Info<< "There are " << nZeroArea
1963  << " faces with area < " << minArea << '.' << nl << endl;
1964  }
1965  else
1966  {
1967  Info<< "All faces have area > " << minArea << '.' << nl << endl;
1968  }
1969  }
1970 
1971  if (nZeroArea > 0)
1972  {
1973  if (report)
1974  {
1976  << nZeroArea << " faces with area < " << minArea
1977  << " found.\n"
1978  << endl;
1979  }
1980 
1981  return true;
1982  }
1983 
1984  return false;
1985 }
1986 
1987 
1990  const bool report,
1991  const scalar warnDet,
1992  const polyMesh& mesh,
1993  const vectorField& faceAreas,
1994  const labelList& checkFaces,
1995  const labelList& affectedCells,
1996  labelHashSet* setPtr
1997 )
1998 {
1999  const cellList& cells = mesh.cells();
2000 
2001  scalar minDet = GREAT;
2002  scalar sumDet = 0.0;
2003  label nSumDet = 0;
2004  label nWarnDet = 0;
2005 
2006  for (const label celli : affectedCells)
2007  {
2008  const cell& cFaces = cells[celli];
2009 
2010  tensor areaSum(Zero);
2011  scalar magAreaSum = 0;
2012 
2013  for (const label facei : cFaces)
2014  {
2015  scalar magArea = mag(faceAreas[facei]);
2016 
2017  magAreaSum += magArea;
2018  areaSum += faceAreas[facei]*(faceAreas[facei]/(magArea+VSMALL));
2019  }
2020 
2021  scalar scaledDet = det(areaSum/(magAreaSum+VSMALL))/0.037037037037037;
2022 
2023  minDet = min(minDet, scaledDet);
2024  sumDet += scaledDet;
2025  ++nSumDet;
2026 
2027  if (scaledDet < warnDet)
2028  {
2029  ++nWarnDet;
2030  if (setPtr)
2031  {
2032  setPtr->insert(cFaces); // Insert all faces of the cell
2033  }
2034  }
2035  }
2036 
2037  reduce(minDet, minOp<scalar>());
2038  reduce(sumDet, sumOp<scalar>());
2039  reduce(nSumDet, sumOp<label>());
2040  reduce(nWarnDet, sumOp<label>());
2041 
2042  if (report)
2043  {
2044  if (nSumDet > 0)
2045  {
2046  Info<< "Cell determinant (1 = uniform cube) : average = "
2047  << sumDet / nSumDet << " min = " << minDet << endl;
2048  }
2049 
2050  if (nWarnDet > 0)
2051  {
2052  Info<< "There are " << nWarnDet
2053  << " cells with determinant < " << warnDet << '.' << nl
2054  << endl;
2055  }
2056  else
2057  {
2058  Info<< "All faces have determinant > " << warnDet << '.' << nl
2059  << endl;
2060  }
2061  }
2062 
2063  if (nWarnDet > 0)
2064  {
2065  if (report)
2066  {
2068  << nWarnDet << " cells with determinant < " << warnDet
2069  << " found.\n"
2070  << endl;
2071  }
2072 
2073  return true;
2074  }
2075 
2076  return false;
2077 }
2078 
2079 
2082  const bool report,
2083  const scalar orthWarn,
2084  const labelList& checkFaces,
2085  const List<labelPair>& baffles,
2086  labelHashSet* setPtr
2087 ) const
2088 {
2089  return checkFaceDotProduct
2090  (
2091  report,
2092  orthWarn,
2093  mesh_,
2094  cellCentres_,
2095  faceAreas_,
2096  checkFaces,
2097  baffles,
2098  setPtr
2099  );
2100 }
2101 
2102 
2105  const bool report,
2106  const scalar minPyrVol,
2107  const pointField& p,
2108  const labelList& checkFaces,
2109  const List<labelPair>& baffles,
2110  labelHashSet* setPtr
2111 ) const
2112 {
2113  return checkFacePyramids
2114  (
2115  report,
2116  minPyrVol,
2117  mesh_,
2118  cellCentres_,
2119  p,
2120  checkFaces,
2121  baffles,
2122  setPtr
2123  );
2124 }
2125 
2126 
2129  const bool report,
2130  const scalar minTetQuality,
2131  const pointField& p,
2132  const labelList& checkFaces,
2133  const List<labelPair>& baffles,
2134  labelHashSet* setPtr
2135 ) const
2136 {
2137  return checkFaceTets
2138  (
2139  report,
2140  minTetQuality,
2141  mesh_,
2142  cellCentres_,
2143  faceCentres_,
2144  p,
2145  checkFaces,
2146  baffles,
2147  setPtr
2148  );
2149 }
2150 
2151 
2152 //bool Foam::polyMeshGeometry::checkFaceSkewness
2153 //(
2154 // const bool report,
2155 // const scalar internalSkew,
2156 // const scalar boundarySkew,
2157 // const labelList& checkFaces,
2158 // const List<labelPair>& baffles,
2159 // labelHashSet* setPtr
2160 //) const
2161 //{
2162 // return checkFaceSkewness
2163 // (
2164 // report,
2165 // internalSkew,
2166 // boundarySkew,
2167 // mesh_,
2168 // cellCentres_,
2169 // faceCentres_,
2170 // faceAreas_,
2171 // checkFaces,
2172 // baffles,
2173 // setPtr
2174 // );
2175 //}
2176 
2177 
2180  const bool report,
2181  const scalar warnWeight,
2182  const labelList& checkFaces,
2183  const List<labelPair>& baffles,
2184  labelHashSet* setPtr
2185 ) const
2186 {
2187  return checkFaceWeights
2188  (
2189  report,
2190  warnWeight,
2191  mesh_,
2192  cellCentres_,
2193  faceCentres_,
2194  faceAreas_,
2195  checkFaces,
2196  baffles,
2197  setPtr
2198  );
2199 }
2200 
2201 
2204  const bool report,
2205  const scalar warnRatio,
2206  const labelList& checkFaces,
2207  const List<labelPair>& baffles,
2208  labelHashSet* setPtr
2209 ) const
2210 {
2211  return checkVolRatio
2212  (
2213  report,
2214  warnRatio,
2215  mesh_,
2216  cellVolumes_,
2217  checkFaces,
2218  baffles,
2219  setPtr
2220  );
2221 }
2222 
2223 
2226  const bool report,
2227  const scalar maxDeg,
2228  const pointField& p,
2229  const labelList& checkFaces,
2230  labelHashSet* setPtr
2231 ) const
2232 {
2233  return checkFaceAngles
2234  (
2235  report,
2236  maxDeg,
2237  mesh_,
2238  faceAreas_,
2239  p,
2240  checkFaces,
2241  setPtr
2242  );
2243 }
2244 
2245 
2248  const bool report,
2249  const scalar minTwist,
2250  const pointField& p,
2251  const labelList& checkFaces,
2252  labelHashSet* setPtr
2253 ) const
2254 {
2255  return checkFaceTwist
2256  (
2257  report,
2258  minTwist,
2259  mesh_,
2260  cellCentres_,
2261  faceAreas_,
2262  faceCentres_,
2263  p,
2264  checkFaces,
2265  setPtr
2266  );
2267 }
2268 
2269 
2272  const bool report,
2273  const scalar minTwist,
2274  const pointField& p,
2275  const labelList& checkFaces,
2276  labelHashSet* setPtr
2277 ) const
2278 {
2279  return checkTriangleTwist
2280  (
2281  report,
2282  minTwist,
2283  mesh_,
2284  faceAreas_,
2285  faceCentres_,
2286  p,
2287  checkFaces,
2288  setPtr
2289  );
2290 }
2291 
2292 
2295  const bool report,
2296  const scalar minFlatness,
2297  const pointField& p,
2298  const labelList& checkFaces,
2299  labelHashSet* setPtr
2300 ) const
2301 {
2302  return checkFaceFlatness
2303  (
2304  report,
2305  minFlatness,
2306  mesh_,
2307  faceAreas_,
2308  faceCentres_,
2309  p,
2310  checkFaces,
2311  setPtr
2312  );
2313 }
2314 
2315 
2318  const bool report,
2319  const scalar minArea,
2320  const labelList& checkFaces,
2321  labelHashSet* setPtr
2322 ) const
2323 {
2324  return checkFaceArea
2325  (
2326  report,
2327  minArea,
2328  mesh_,
2329  faceAreas_,
2330  checkFaces,
2331  setPtr
2332  );
2333 }
2334 
2335 
2338  const bool report,
2339  const scalar warnDet,
2340  const labelList& checkFaces,
2341  const labelList& affectedCells,
2342  labelHashSet* setPtr
2343 ) const
2344 {
2345  return checkCellDeterminant
2346  (
2347  report,
2348  warnDet,
2349  mesh_,
2350  faceAreas_,
2351  checkFaces,
2352  affectedCells,
2353  setPtr
2354  );
2355 }
2356 
2357 
2358 // ************************************************************************* //
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::polyMeshGeometry::checkFaceDotProduct
static bool checkFaceDotProduct(const bool report, const scalar orthWarn, const polyMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
See primitiveMesh.
Definition: polyMeshGeometry.C:360
Foam::Tensor< scalar >
Foam::maxOp
Definition: ops.H:223
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::polyBoundaryMesh
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Definition: polyBoundaryMesh.H:63
primitiveMeshTools.H
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
polyMeshGeometry.H
Foam::minOp
Definition: ops.H:224
Foam::sin
dimensionedScalar sin(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:264
Foam::boundBox::invertedBox
static const boundBox invertedBox
A large inverted boundBox: min/max == +/- ROOTVGREAT.
Definition: boundBox.H:86
Foam::polyMeshGeometry::checkFaceWeights
static bool checkFaceWeights(const bool report, const scalar warnWeight, const polyMesh &mesh, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
Interpolation weights (0.5 for regular mesh)
Definition: polyMeshGeometry.C:1135
Foam::primitiveMesh::cells
const cellList & cells() const
Definition: primitiveMeshCells.C:138
unitConversion.H
Unit conversion functions.
Foam::syncTools::swapBoundaryFaceList
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &faceValues)
Swap coupled boundary face values. Uses eqOp.
Definition: syncTools.H:445
tetPointRef.H
Foam::polyMeshGeometry::correct
void correct()
Take over properties from mesh.
Definition: polyMeshGeometry.C:337
Foam::polyMesh::boundaryMesh
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:444
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::polyMeshGeometry::checkFaceTets
static bool checkFaceTets(const bool report, const scalar minPyrVol, const polyMesh &, const vectorField &cellCentres, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *)
See primitiveMesh.
Definition: polyMeshGeometry.C:709
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::polyMeshGeometry::polyMeshGeometry
polyMeshGeometry(const polyMesh &)
Construct from mesh.
Definition: polyMeshGeometry.C:327
Foam::HashSet< label, Hash< label > >
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
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
polyMeshTetDecomposition.H
Foam::sumOp
Definition: ops.H:213
Foam::tetPointRef
tetrahedron< point, const point & > tetPointRef
A tetrahedron using referred points.
Definition: tetPointRef.H:47
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::polyMeshGeometry::checkFaceSkewness
static bool checkFaceSkewness(const bool report, const scalar internalSkew, const scalar boundarySkew, const polyMesh &mesh, const pointField &points, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
See primitiveMesh.
Definition: polyMeshGeometry.C:930
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
Foam::polyMeshTetDecomposition::findSharedBasePoint
static label findSharedBasePoint(const polyMesh &mesh, label fI, const point &nCc, scalar tol, bool report=false)
Find the first suitable base point to use for a minimum.
Definition: polyMeshTetDecomposition.C:134
Foam::triangle
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:59
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:51
Foam::polyMeshTetDecomposition::findBasePoint
static label findBasePoint(const polyMesh &mesh, label fI, scalar tol, bool report=false)
Find the base point to use for a minimum triangle.
Definition: polyMeshTetDecomposition.C:189
Foam::polyMeshGeometry::affectedCells
static labelList affectedCells(const polyMesh &, const labelList &changedFaces)
Helper function: get affected cells from faces.
Definition: polyMeshGeometry.C:183
Foam::polyMeshGeometry::checkFaceArea
static bool checkFaceArea(const bool report, const scalar minArea, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
Small faces.
Definition: polyMeshGeometry.C:1932
Foam::Field< vector >
Foam::syncTools::swapBoundaryCellPositions
static void swapBoundaryCellPositions(const polyMesh &mesh, const UList< point > &cellData, List< point > &neighbourCellData)
Swap to obtain neighbour cell positions for all boundary faces.
Definition: syncTools.C:34
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::polyMesh::faceOwner
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1107
SeriousErrorInFunction
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
Definition: messageStream.H:306
Foam::cellList
List< cell > cellList
A List of cells.
Definition: cellListFwd.H:47
Foam::primitiveMesh::nBoundaryFaces
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
Definition: primitiveMeshI.H:84
Foam::polyBoundaryMesh::whichPatch
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Definition: polyBoundaryMesh.C:812
Foam::radToDeg
constexpr scalar radToDeg(const scalar rad) noexcept
Conversion from radians to degrees.
Definition: unitConversion.H:54
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
Foam::syncTools::swapBoundaryFacePositions
static void swapBoundaryFacePositions(const polyMesh &mesh, UList< point > &positions)
Swap coupled positions. Uses eqOp.
Definition: syncTools.H:461
Foam::FatalError
error FatalError
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::degToRad
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Definition: unitConversion.H:48
Foam::polyMeshGeometry::checkFaceTwist
static bool checkFaceTwist(const bool report, const scalar minTwist, const polyMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Triangle (from face-centre decomposition) normal v.s.
Definition: polyMeshGeometry.C:1547
Foam::normalised
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:487
Foam::polyMeshGeometry::checkFaceFlatness
static bool checkFaceFlatness(const bool report, const scalar minFlatness, const polyMesh &, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Area of faces v.s. sum of triangle areas.
Definition: polyMeshGeometry.C:1839
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::primitiveMesh::isInternalFace
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
Definition: primitiveMeshI.H:103
Foam::polyMesh::faces
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1094
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::Pair< label >
f
labelList f(nPoints)
Foam::faceList
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
Foam::Vector< scalar >
Foam::primitiveMesh::nInternalFaces
label nInternalFaces() const noexcept
Number of internal faces.
Definition: primitiveMeshI.H:78
Foam::List< label >
Foam::acos
dimensionedScalar acos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:268
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::HashSet::insert
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
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
patches
const polyBoundaryMesh & patches
Definition: convertProcessorPatches.H:65
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::VectorSpace< Vector< scalar >, scalar, 3 >::zero
static const Vector< scalar > zero
Definition: VectorSpace.H:115
Foam::primitiveMesh::nFaces
label nFaces() const noexcept
Number of mesh faces.
Definition: primitiveMeshI.H:90
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:72
cells
const cellShapeList & cells
Definition: gmvOutputHeader.H:3
Foam::det
dimensionedScalar det(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:62
Foam::polyMeshGeometry::checkFacePyramids
static bool checkFacePyramids(const bool report, const scalar minPyrVol, const polyMesh &, const vectorField &cellCentres, const pointField &p, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *)
See primitiveMesh.
Definition: polyMeshGeometry.C:538
coupled
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
Foam::VectorSpace::size
static constexpr direction size() noexcept
The number of elements in the VectorSpace = Ncmpts.
Definition: VectorSpace.H:176
Foam::point
vector point
Point is a vector.
Definition: point.H:43
Foam::polyMeshGeometry::checkCellDeterminant
static bool checkCellDeterminant(const bool report, const scalar minDet, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, const labelList &affectedCells, labelHashSet *setPtr)
Area of internal faces v.s. boundary faces.
Definition: polyMeshGeometry.C:1989
Foam::labelHashSet
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
Foam::polyMeshGeometry::checkTriangleTwist
static bool checkTriangleTwist(const bool report, const scalar minTwist, const polyMesh &, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Consecutive triangle (from face-centre decomposition) normals.
Definition: polyMeshGeometry.C:1693
Foam::pyramidPointFaceRef
pyramid< point, const point &, const face & > pyramidPointFaceRef
A pyramid using referred points and faces.
Definition: pyramidPointFaceRef.H:48
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
pyramidPointFaceRef.H
Foam::cell
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:54
Foam::asin
dimensionedScalar asin(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:267
Foam::polyMesh::faceNeighbour
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1113
Foam::polyMeshGeometry::checkVolRatio
static bool checkVolRatio(const bool report, const scalar warnRatio, const polyMesh &mesh, const scalarField &cellVolumes, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
Cell volume ratio of neighbouring cells (1 for regular mesh)
Definition: polyMeshGeometry.C:1283
Foam::triPointRef
triangle< point, const point & > triPointRef
A triangle using referred points.
Definition: triangle.H:71
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:265
Foam::polyMeshGeometry::checkFaceAngles
static bool checkFaceAngles(const bool report, const scalar maxDeg, const polyMesh &mesh, const vectorField &faceAreas, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
See primitiveMesh.
Definition: polyMeshGeometry.C:1419