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