primitiveMeshCheck.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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 "primitiveMesh.H"
27 #include "pyramidPointFaceRef.H"
28 #include "ListOps.H"
29 #include "unitConversion.H"
30 #include "SortableList.H"
31 #include "EdgeMap.H"
32 #include "primitiveMeshTools.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 Foam::scalar Foam::primitiveMesh::closedThreshold_ = 1.0e-6;
37 Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000;
38 Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70; // deg
40 Foam::scalar Foam::primitiveMesh::planarCosAngle_ = 1.0e-6;
41 
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
46 (
47  const vectorField& areas,
48  const bool report,
49  const PackedBoolList& internalOrCoupledFaces
50 ) const
51 {
52  if (debug)
53  {
55  << "Checking whether the boundary is closed" << endl;
56  }
57 
58  // Loop through all boundary faces and sum up the face area vectors.
59  // For a closed boundary, this should be zero in all vector components
60 
61  vector sumClosed(Zero);
62  scalar sumMagClosedBoundary = 0;
63 
64  for (label facei = nInternalFaces(); facei < areas.size(); facei++)
65  {
66  if (!internalOrCoupledFaces.size() || !internalOrCoupledFaces[facei])
67  {
68  sumClosed += areas[facei];
69  sumMagClosedBoundary += mag(areas[facei]);
70  }
71  }
72 
73  reduce(sumClosed, sumOp<vector>());
74  reduce(sumMagClosedBoundary, sumOp<scalar>());
75 
76  vector openness = sumClosed/(sumMagClosedBoundary + vSmall);
77 
78  if (cmptMax(cmptMag(openness)) > closedThreshold_)
79  {
80  if (debug || report)
81  {
82  Info<< " ***Boundary openness " << openness
83  << " possible hole in boundary description."
84  << endl;
85  }
86 
87  return true;
88  }
89  else
90  {
91  if (debug || report)
92  {
93  Info<< " Boundary openness " << openness << " OK."
94  << endl;
95  }
96 
97  return false;
98  }
99 }
100 
101 
103 (
104  const vectorField& faceAreas,
105  const scalarField& cellVolumes,
106  const bool report,
107  labelHashSet* setPtr,
108  labelHashSet* aspectSetPtr,
109  const Vector<label>& meshD
110 ) const
111 {
112  if (debug)
113  {
115  << "Checking whether cells are closed" << endl;
116  }
117 
118  // Check that all cells labels are valid
119  const cellList& c = cells();
120 
121  label nErrorClosed = 0;
122 
123  forAll(c, cI)
124  {
125  const cell& curCell = c[cI];
126 
127  if (min(curCell) < 0 || max(curCell) > nFaces())
128  {
129  if (setPtr)
130  {
131  setPtr->insert(cI);
132  }
133 
134  nErrorClosed++;
135  }
136  }
137 
138  if (nErrorClosed > 0)
139  {
140  if (debug || report)
141  {
142  Info<< " ***Cells with invalid face labels found, number of cells "
143  << nErrorClosed << endl;
144  }
145 
146  return true;
147  }
148 
149 
150  scalarField openness;
151  scalarField aspectRatio;
152  primitiveMeshTools::cellClosedness
153  (
154  *this,
155  meshD,
156  faceAreas,
157  cellVolumes,
158  openness,
159  aspectRatio
160  );
161 
162  label nOpen = 0;
163  scalar maxOpennessCell = max(openness);
164  label nAspect = 0;
165  scalar maxAspectRatio = max(aspectRatio);
166 
167  // Check the sums
168  forAll(openness, celli)
169  {
170  if (openness[celli] > closedThreshold_)
171  {
172  if (setPtr)
173  {
174  setPtr->insert(celli);
175  }
176 
177  nOpen++;
178  }
179 
180  if (aspectRatio[celli] > aspectThreshold_)
181  {
182  if (aspectSetPtr)
183  {
184  aspectSetPtr->insert(celli);
185  }
186 
187  nAspect++;
188  }
189  }
190 
191  reduce(nOpen, sumOp<label>());
192  reduce(maxOpennessCell, maxOp<scalar>());
193 
194  reduce(nAspect, sumOp<label>());
195  reduce(maxAspectRatio, maxOp<scalar>());
196 
197 
198  if (nOpen > 0)
199  {
200  if (debug || report)
201  {
202  Info<< " ***Open cells found, max cell openness: "
203  << maxOpennessCell << ", number of open cells " << nOpen
204  << endl;
205  }
206 
207  return true;
208  }
209 
210  if (nAspect > 0)
211  {
212  if (debug || report)
213  {
214  Info<< " ***High aspect ratio cells found, Max aspect ratio: "
215  << maxAspectRatio
216  << ", number of cells " << nAspect
217  << endl;
218  }
219 
220  return true;
221  }
222 
223  if (debug || report)
224  {
225  Info<< " Max cell openness = " << maxOpennessCell << " OK." << nl
226  << " Max aspect ratio = " << maxAspectRatio << " OK."
227  << endl;
228  }
229 
230  return false;
231 }
232 
233 
235 (
236  const vectorField& faceAreas,
237  const bool report,
238  const bool detailedReport,
239  labelHashSet* setPtr
240 ) const
241 {
242  if (debug)
243  {
244  InfoInFunction << "Checking face area magnitudes" << endl;
245  }
246 
247  const scalarField magFaceAreas(mag(faceAreas));
248 
249  scalar minArea = great;
250  scalar maxArea = -great;
251 
252  forAll(magFaceAreas, facei)
253  {
254  if (magFaceAreas[facei] < vSmall)
255  {
256  if (setPtr)
257  {
258  setPtr->insert(facei);
259  }
260  if (detailedReport)
261  {
262  if (isInternalFace(facei))
263  {
264  Pout<< "Zero or negative face area detected for "
265  << "internal face "<< facei << " between cells "
266  << faceOwner()[facei] << " and "
267  << faceNeighbour()[facei]
268  << ". Face area magnitude = " << magFaceAreas[facei]
269  << endl;
270  }
271  else
272  {
273  Pout<< "Zero or negative face area detected for "
274  << "boundary face " << facei << " next to cell "
275  << faceOwner()[facei] << ". Face area magnitude = "
276  << magFaceAreas[facei] << endl;
277  }
278  }
279  }
280 
281  minArea = min(minArea, magFaceAreas[facei]);
282  maxArea = max(maxArea, magFaceAreas[facei]);
283  }
284 
285  reduce(minArea, minOp<scalar>());
286  reduce(maxArea, maxOp<scalar>());
287 
288  if (minArea < vSmall)
289  {
290  if (debug || report)
291  {
292  Info<< " ***Zero or negative face area detected. "
293  "Minimum area: " << minArea << endl;
294  }
295 
296  return true;
297  }
298  else
299  {
300  if (debug || report)
301  {
302  Info<< " Minimum face area = " << minArea
303  << ". Maximum face area = " << maxArea
304  << ". Face area magnitudes OK." << endl;
305  }
306 
307  return false;
308  }
309 }
310 
311 
313 (
314  const scalarField& vols,
315  const bool report,
316  const bool detailedReport,
317  labelHashSet* setPtr
318 ) const
319 {
320  if (debug)
321  {
322  InfoInFunction << "Checking cell volumes" << endl;
323  }
324 
325  scalar minVolume = great;
326  scalar maxVolume = -great;
327 
328  label nNegVolCells = 0;
329 
330  forAll(vols, celli)
331  {
332  if (vols[celli] < vSmall)
333  {
334  if (setPtr)
335  {
336  setPtr->insert(celli);
337  }
338  if (detailedReport)
339  {
340  Pout<< "Zero or negative cell volume detected for cell "
341  << celli << ". Volume = " << vols[celli] << endl;
342  }
343 
344  nNegVolCells++;
345  }
346 
347  minVolume = min(minVolume, vols[celli]);
348  maxVolume = max(maxVolume, vols[celli]);
349  }
350 
351  reduce(minVolume, minOp<scalar>());
352  reduce(maxVolume, maxOp<scalar>());
353  reduce(nNegVolCells, sumOp<label>());
354 
355  if (minVolume < vSmall)
356  {
357  if (debug || report)
358  {
359  Info<< " ***Zero or negative cell volume detected. "
360  << "Minimum negative volume: " << minVolume
361  << ", Number of negative volume cells: " << nNegVolCells
362  << endl;
363  }
364 
365  return true;
366  }
367  else
368  {
369  if (debug || report)
370  {
371  Info<< " Min volume = " << minVolume
372  << ". Max volume = " << maxVolume
373  << ". Total volume = " << gSum(vols)
374  << ". Cell volumes OK." << endl;
375  }
376 
377  return false;
378  }
379 }
380 
381 
383 (
384  const vectorField& fAreas,
385  const vectorField& cellCtrs,
386  const bool report,
387  labelHashSet* setPtr
388 ) const
389 {
390  if (debug)
391  {
392  InfoInFunction << "Checking mesh non-orthogonality" << endl;
393  }
394 
395 
396  tmp<scalarField> tortho = primitiveMeshTools::faceOrthogonality
397  (
398  *this,
399  fAreas,
400  cellCtrs
401  );
402  const scalarField& ortho = tortho();
403 
404  // Severe nonorthogonality threshold
405  const scalar severeNonorthogonalityThreshold =
406  ::cos(degToRad(nonOrthThreshold_));
407 
408  scalar minDDotS = min(ortho);
409 
410  scalar sumDDotS = sum(ortho);
411 
412  label severeNonOrth = 0;
413 
414  label errorNonOrth = 0;
415 
416 
417  forAll(ortho, facei)
418  {
419  if (ortho[facei] < severeNonorthogonalityThreshold)
420  {
421  if (ortho[facei] > small)
422  {
423  if (setPtr)
424  {
425  setPtr->insert(facei);
426  }
427 
428  severeNonOrth++;
429  }
430  else
431  {
432  if (setPtr)
433  {
434  setPtr->insert(facei);
435  }
436 
437  errorNonOrth++;
438  }
439  }
440  }
441 
442  reduce(minDDotS, minOp<scalar>());
443  reduce(sumDDotS, sumOp<scalar>());
444  reduce(severeNonOrth, sumOp<label>());
445  reduce(errorNonOrth, sumOp<label>());
446 
447  if (debug || report)
448  {
449  label neiSize = ortho.size();
450  reduce(neiSize, sumOp<label>());
451 
452  if (neiSize > 0)
453  {
454  if (debug || report)
455  {
456  Info<< " Mesh non-orthogonality Max: "
457  << radToDeg(::acos(minDDotS))
458  << " average: " << radToDeg(::acos(sumDDotS/neiSize))
459  << endl;
460  }
461  }
462 
463  if (severeNonOrth > 0)
464  {
465  Info<< " *Number of severely non-orthogonal faces: "
466  << severeNonOrth << "." << endl;
467  }
468  }
469 
470  if (errorNonOrth > 0)
471  {
472  if (debug || report)
473  {
474  Info<< " ***Number of non-orthogonality errors: "
475  << errorNonOrth << "." << endl;
476  }
477 
478  return true;
479  }
480  else
481  {
482  if (debug || report)
483  {
484  Info<< " Non-orthogonality check OK." << endl;
485  }
486 
487  return false;
488  }
489 }
490 
491 
493 (
494  const pointField& points,
495  const vectorField& ctrs,
496  const bool report,
497  const bool detailedReport,
498  const scalar minPyrVol,
499  labelHashSet* setPtr
500 ) const
501 {
502  if (debug)
503  {
504  InfoInFunction << "Checking face orientation" << endl;
505  }
506 
507  const labelList& own = faceOwner();
508  const labelList& nei = faceNeighbour();
509  const faceList& f = faces();
510 
511 
512  scalarField ownPyrVol;
513  scalarField neiPyrVol;
514  primitiveMeshTools::facePyramidVolume
515  (
516  *this,
517  points,
518  ctrs,
519  ownPyrVol,
520  neiPyrVol
521  );
522 
523 
524  label nErrorPyrs = 0;
525 
526  forAll(ownPyrVol, facei)
527  {
528  if (ownPyrVol[facei] < minPyrVol)
529  {
530  if (setPtr)
531  {
532  setPtr->insert(facei);
533  }
534  if (detailedReport)
535  {
536  Pout<< "Negative pyramid volume: " << ownPyrVol[facei]
537  << " for face " << facei << " " << f[facei]
538  << " and owner cell: " << own[facei] << endl
539  << "Owner cell vertex labels: "
540  << cells()[own[facei]].labels(faces())
541  << endl;
542  }
543 
544  nErrorPyrs++;
545  }
546 
547  if (isInternalFace(facei))
548  {
549  if (neiPyrVol[facei] < minPyrVol)
550  {
551  if (setPtr)
552  {
553  setPtr->insert(facei);
554  }
555  if (detailedReport)
556  {
557  Pout<< "Negative pyramid volume: " << neiPyrVol[facei]
558  << " for face " << facei << " " << f[facei]
559  << " and neighbour cell: " << nei[facei] << nl
560  << "Neighbour cell vertex labels: "
561  << cells()[nei[facei]].labels(faces())
562  << endl;
563  }
564  nErrorPyrs++;
565  }
566  }
567  }
568 
569  reduce(nErrorPyrs, sumOp<label>());
570 
571  if (nErrorPyrs > 0)
572  {
573  if (debug || report)
574  {
575  Info<< " ***Error in face pyramids: "
576  << nErrorPyrs << " faces are incorrectly oriented."
577  << endl;
578  }
579 
580  return true;
581  }
582  else
583  {
584  if (debug || report)
585  {
586  Info<< " Face pyramids OK." << endl;
587  }
588 
589  return false;
590  }
591 }
592 
593 
595 (
596  const pointField& points,
597  const vectorField& fCtrs,
598  const vectorField& fAreas,
599  const vectorField& cellCtrs,
600  const bool report,
601  labelHashSet* setPtr
602 ) const
603 {
604  if (debug)
605  {
606  InfoInFunction << "Checking face skewness" << endl;
607  }
608 
609  // Warn if the skew correction vector is more than skewWarning times
610  // larger than the face area vector
611 
612  tmp<scalarField> tskewness = primitiveMeshTools::faceSkewness
613  (
614  *this,
615  points,
616  fCtrs,
617  fAreas,
618  cellCtrs
619  );
620  const scalarField& skewness = tskewness();
621 
622  scalar maxSkew = max(skewness);
623  label nWarnSkew = 0;
624 
625  forAll(skewness, facei)
626  {
627  // Check if the skewness vector is greater than the PN vector.
628  // This does not cause trouble but is a good indication of a poor mesh.
629  if (skewness[facei] > skewThreshold_)
630  {
631  if (setPtr)
632  {
633  setPtr->insert(facei);
634  }
635 
636  nWarnSkew++;
637  }
638  }
639 
640  reduce(maxSkew, maxOp<scalar>());
641  reduce(nWarnSkew, sumOp<label>());
642 
643  if (nWarnSkew > 0)
644  {
645  if (debug || report)
646  {
647  Info<< " ***Max skewness = " << maxSkew
648  << ", " << nWarnSkew << " highly skew faces detected"
649  " which may impair the quality of the results"
650  << endl;
651  }
652 
653  return true;
654  }
655  else
656  {
657  if (debug || report)
658  {
659  Info<< " Max skewness = " << maxSkew << " OK." << endl;
660  }
661 
662  return false;
663  }
664 }
665 
666 
668 (
669  const pointField& points,
670  const vectorField& faceAreas,
671  const bool report,
672  const scalar maxDeg,
673  labelHashSet* setPtr
674 ) const
675 {
676  if (debug)
677  {
678  InfoInFunction << "Checking face angles" << endl;
679  }
680 
681  if (maxDeg < -small || maxDeg > 180+small)
682  {
684  << "maxDeg should be [0..180] but is now " << maxDeg
685  << exit(FatalError);
686  }
687 
688  const scalar maxSin = Foam::sin(degToRad(maxDeg));
689 
690 
691  tmp<scalarField> tfaceAngles = primitiveMeshTools::faceConcavity
692  (
693  maxSin,
694  *this,
695  points,
696  faceAreas
697  );
698  const scalarField& faceAngles = tfaceAngles();
699 
700  scalar maxEdgeSin = max(faceAngles);
701 
702  label nConcave = 0;
703 
704  forAll(faceAngles, facei)
705  {
706  if (faceAngles[facei] > small)
707  {
708  nConcave++;
709 
710  if (setPtr)
711  {
712  setPtr->insert(facei);
713  }
714  }
715  }
716 
717  reduce(nConcave, sumOp<label>());
718  reduce(maxEdgeSin, maxOp<scalar>());
719 
720  if (nConcave > 0)
721  {
722  scalar maxConcaveDegr =
723  radToDeg(Foam::asin(Foam::min(1.0, maxEdgeSin)));
724 
725  if (debug || report)
726  {
727  Info<< " *There are " << nConcave
728  << " faces with concave angles between consecutive"
729  << " edges. Max concave angle = " << maxConcaveDegr
730  << " degrees." << endl;
731  }
732 
733  return true;
734  }
735  else
736  {
737  if (debug || report)
738  {
739  Info<< " All angles in faces OK." << endl;
740  }
741 
742  return false;
743  }
744 }
745 
746 
748 (
749  const pointField& points,
750  const vectorField& faceCentres,
751  const vectorField& faceAreas,
752  const bool report,
753  const scalar warnFlatness,
754  labelHashSet* setPtr
755 ) const
756 {
757  if (debug)
758  {
759  InfoInFunction << "Checking face flatness" << endl;
760  }
761 
762  if (warnFlatness < 0 || warnFlatness > 1)
763  {
765  << "warnFlatness should be [0..1] but is now " << warnFlatness
766  << exit(FatalError);
767  }
768 
769  const faceList& fcs = faces();
770 
771  tmp<scalarField> tfaceFlatness = primitiveMeshTools::faceFlatness
772  (
773  *this,
774  points,
775  faceCentres,
776  faceAreas
777  );
778  const scalarField& faceFlatness = tfaceFlatness();
779 
780  scalarField magAreas(mag(faceAreas));
781 
782  scalar minFlatness = great;
783  scalar sumFlatness = 0;
784  label nSummed = 0;
785  label nWarped = 0;
786 
787  forAll(faceFlatness, facei)
788  {
789  if (fcs[facei].size() > 3 && magAreas[facei] > vSmall)
790  {
791  sumFlatness += faceFlatness[facei];
792  nSummed++;
793 
794  minFlatness = min(minFlatness, faceFlatness[facei]);
795 
796  if (faceFlatness[facei] < warnFlatness)
797  {
798  nWarped++;
799 
800  if (setPtr)
801  {
802  setPtr->insert(facei);
803  }
804  }
805  }
806  }
807 
808 
809  reduce(nWarped, sumOp<label>());
810  reduce(minFlatness, minOp<scalar>());
811 
812  reduce(nSummed, sumOp<label>());
813  reduce(sumFlatness, sumOp<scalar>());
814 
815  if (debug || report)
816  {
817  if (nSummed > 0)
818  {
819  Info<< " Face flatness (1 = flat, 0 = butterfly) : min = "
820  << minFlatness << " average = " << sumFlatness / nSummed
821  << endl;
822  }
823  }
824 
825 
826  if (nWarped> 0)
827  {
828  if (debug || report)
829  {
830  Info<< " *There are " << nWarped
831  << " faces with ratio between projected and actual area < "
832  << warnFlatness << endl;
833 
834  Info<< " Minimum ratio (minimum flatness, maximum warpage) = "
835  << minFlatness << endl;
836  }
837 
838  return true;
839  }
840  else
841  {
842  if (debug || report)
843  {
844  Info<< " All face flatness OK." << endl;
845  }
846 
847  return false;
848  }
849 }
850 
851 
853 (
854  const vectorField& fAreas,
855  const pointField& fCentres,
856  const bool report,
857  labelHashSet* setPtr
858 ) const
859 {
860  if (debug)
861  {
862  InfoInFunction << "Checking for concave cells" << endl;
863  }
864 
865  const cellList& c = cells();
866  const labelList& fOwner = faceOwner();
867 
868  label nConcaveCells = 0;
869 
870  forAll(c, celli)
871  {
872  const cell& cFaces = c[celli];
873 
874  bool concave = false;
875 
876  forAll(cFaces, i)
877  {
878  if (concave)
879  {
880  break;
881  }
882 
883  label fI = cFaces[i];
884 
885  const point& fC = fCentres[fI];
886 
887  vector fN = fAreas[fI];
888 
889  fN /= max(mag(fN), vSmall);
890 
891  // Flip normal if required so that it is always pointing out of
892  // the cell
893  if (fOwner[fI] != celli)
894  {
895  fN *= -1;
896  }
897 
898  // Is the centre of any other face of the cell on the
899  // wrong side of the plane of this face?
900 
901  forAll(cFaces, j)
902  {
903  if (j != i)
904  {
905  label fJ = cFaces[j];
906 
907  const point& pt = fCentres[fJ];
908 
909  // If the cell is concave, the point will be on the
910  // positive normal side of the plane of f, defined by
911  // its centre and normal, and the angle between (pt -
912  // fC) and fN will be less than 90 degrees, so the dot
913  // product will be positive.
914 
915  vector pC = (pt - fC);
916 
917  pC /= max(mag(pC), vSmall);
918 
919  if ((pC & fN) > -planarCosAngle_)
920  {
921  // Concave or planar face
922 
923  concave = true;
924 
925  if (setPtr)
926  {
927  setPtr->insert(celli);
928  }
929 
930  nConcaveCells++;
931 
932  break;
933  }
934  }
935  }
936  }
937  }
938 
939  reduce(nConcaveCells, sumOp<label>());
940 
941  if (nConcaveCells > 0)
942  {
943  if (debug || report)
944  {
945  Info<< " ***Concave cells (using face planes) found,"
946  << " number of cells: " << nConcaveCells << endl;
947  }
948 
949  return true;
950  }
951  else
952  {
953  if (debug || report)
954  {
955  Info<< " Concave cell check OK." << endl;
956  }
957 
958  return false;
959  }
960 
961  return false;
962 }
963 
964 
966 (
967  const bool report,
968  labelHashSet* setPtr
969 ) const
970 {
971  if (debug)
972  {
973  InfoInFunction << "Checking face ordering" << endl;
974  }
975 
976  // Check whether internal faces are ordered in the upper triangular order
977  const labelList& own = faceOwner();
978  const labelList& nei = faceNeighbour();
979 
980  const cellList& c = cells();
981 
982  label internal = nInternalFaces();
983 
984  // Has error occurred?
985  bool error = false;
986  // Have multiple faces been detected?
987  label nMultipleCells = false;
988 
989  // Loop through faceCells once more and make sure that for internal cell
990  // the first label is smaller
991  for (label facei = 0; facei < internal; facei++)
992  {
993  if (own[facei] >= nei[facei])
994  {
995  error = true;
996 
997  if (setPtr)
998  {
999  setPtr->insert(facei);
1000  }
1001  }
1002  }
1003 
1004  // Loop through all cells. For each cell, find the face that is internal
1005  // and add it to the check list (upper triangular order).
1006  // Once the list is completed, check it against the faceCell list
1007 
1008  forAll(c, celli)
1009  {
1010  const labelList& curFaces = c[celli];
1011 
1012  // Neighbouring cells
1013  SortableList<label> nbr(curFaces.size());
1014 
1015  forAll(curFaces, i)
1016  {
1017  label facei = curFaces[i];
1018 
1019  if (facei >= nInternalFaces())
1020  {
1021  // Sort last
1022  nbr[i] = labelMax;
1023  }
1024  else
1025  {
1026  label nbrCelli = nei[facei];
1027 
1028  if (nbrCelli == celli)
1029  {
1030  nbrCelli = own[facei];
1031  }
1032 
1033  if (celli < nbrCelli)
1034  {
1035  // celli is master
1036  nbr[i] = nbrCelli;
1037  }
1038  else
1039  {
1040  // nbrCell is master. Let it handle this face.
1041  nbr[i] = labelMax;
1042  }
1043  }
1044  }
1045 
1046  nbr.sort();
1047 
1048  // Now nbr holds the cellCells in incremental order. Check:
1049  // - neighbouring cells appear only once. Since nbr is sorted this
1050  // is simple check on consecutive elements
1051  // - faces indexed in same order as nbr are incrementing as well.
1052 
1053  label prevCell = nbr[0];
1054  label prevFace = curFaces[nbr.indices()[0]];
1055 
1056  bool hasMultipleFaces = false;
1057 
1058  for (label i = 1; i < nbr.size(); i++)
1059  {
1060  label thisCell = nbr[i];
1061  label thisFace = curFaces[nbr.indices()[i]];
1062 
1063  if (thisCell == labelMax)
1064  {
1065  break;
1066  }
1067 
1068  if (thisCell == prevCell)
1069  {
1070  hasMultipleFaces = true;
1071 
1072  if (setPtr)
1073  {
1074  setPtr->insert(prevFace);
1075  setPtr->insert(thisFace);
1076  }
1077  }
1078  else if (thisFace < prevFace)
1079  {
1080  error = true;
1081 
1082  if (setPtr)
1083  {
1084  setPtr->insert(thisFace);
1085  }
1086  }
1087 
1088  prevCell = thisCell;
1089  prevFace = thisFace;
1090  }
1091 
1092  if (hasMultipleFaces)
1093  {
1094  nMultipleCells++;
1095  }
1096  }
1097 
1098  reduce(error, orOp<bool>());
1099  reduce(nMultipleCells, sumOp<label>());
1100 
1101  if ((debug || report) && nMultipleCells > 0)
1102  {
1103  Info<< " <<Found " << nMultipleCells
1104  << " neighbouring cells with multiple in between faces." << endl;
1105  }
1106 
1107  if (error)
1108  {
1109  if (debug || report)
1110  {
1111  Info<< " ***Faces not in upper triangular order." << endl;
1112  }
1113 
1114  return true;
1115  }
1116  else
1117  {
1118  if (debug || report)
1119  {
1120  Info<< " Upper triangular ordering OK." << endl;
1121  }
1122 
1123  return false;
1124  }
1125 }
1126 
1127 
1130  const bool report,
1131  labelHashSet* setPtr
1132 ) const
1133 {
1134  if (debug)
1135  {
1136  InfoInFunction << "Checking topological cell openness" << endl;
1137  }
1138 
1139  label nOpenCells = 0;
1140 
1141  const faceList& f = faces();
1142  const cellList& c = cells();
1143 
1144  forAll(c, celli)
1145  {
1146  const labelList& curFaces = c[celli];
1147 
1148  const edgeList cellEdges = c[celli].edges(f);
1149 
1150  labelList edgeUsage(cellEdges.size(), 0);
1151 
1152  forAll(curFaces, facei)
1153  {
1154  edgeList curFaceEdges = f[curFaces[facei]].edges();
1155 
1156  forAll(curFaceEdges, faceEdgeI)
1157  {
1158  const edge& curEdge = curFaceEdges[faceEdgeI];
1159 
1160  forAll(cellEdges, cellEdgeI)
1161  {
1162  if (cellEdges[cellEdgeI] == curEdge)
1163  {
1164  edgeUsage[cellEdgeI]++;
1165  break;
1166  }
1167  }
1168  }
1169  }
1170 
1171  edgeList singleEdges(cellEdges.size());
1172  label nSingleEdges = 0;
1173 
1174  forAll(edgeUsage, edgeI)
1175  {
1176  if (edgeUsage[edgeI] == 1)
1177  {
1178  singleEdges[nSingleEdges] = cellEdges[edgeI];
1179  nSingleEdges++;
1180  }
1181  else if (edgeUsage[edgeI] != 2)
1182  {
1183  if (setPtr)
1184  {
1185  setPtr->insert(celli);
1186  }
1187  }
1188  }
1189 
1190  if (nSingleEdges > 0)
1191  {
1192  if (setPtr)
1193  {
1194  setPtr->insert(celli);
1195  }
1196 
1197  nOpenCells++;
1198  }
1199  }
1200 
1201  reduce(nOpenCells, sumOp<label>());
1202 
1203  if (nOpenCells > 0)
1204  {
1205  if (debug || report)
1206  {
1207  Info<< " ***Open cells found, number of cells: " << nOpenCells
1208  << ". This problem may be fixable using the zipUpMesh utility."
1209  << endl;
1210  }
1211 
1212  return true;
1213  }
1214  else
1215  {
1216  if (debug || report)
1217  {
1218  Info<< " Topological cell zip-up check OK." << endl;
1219  }
1220 
1221  return false;
1222  }
1223 }
1224 
1225 
1228  const bool report,
1229  labelHashSet* setPtr
1230 ) const
1231 {
1232  if (debug)
1233  {
1234  InfoInFunction << "Checking face vertices" << endl;
1235  }
1236 
1237  // Check that all vertex labels are valid
1238  const faceList& f = faces();
1239 
1240  label nErrorFaces = 0;
1241 
1242  forAll(f, fI)
1243  {
1244  const face& curFace = f[fI];
1245 
1246  if (min(curFace) < 0 || max(curFace) > nPoints())
1247  {
1248  if (setPtr)
1249  {
1250  setPtr->insert(fI);
1251  }
1252 
1253  nErrorFaces++;
1254  }
1255 
1256  // Uniqueness of vertices
1257  labelHashSet facePoints(2*curFace.size());
1258 
1259  forAll(curFace, fp)
1260  {
1261  bool inserted = facePoints.insert(curFace[fp]);
1262 
1263  if (!inserted)
1264  {
1265  if (setPtr)
1266  {
1267  setPtr->insert(fI);
1268  }
1269 
1270  nErrorFaces++;
1271  }
1272  }
1273  }
1274 
1275  reduce(nErrorFaces, sumOp<label>());
1276 
1277  if (nErrorFaces > 0)
1278  {
1279  if (debug || report)
1280  {
1281  Info<< " Faces with invalid vertex labels found, "
1282  << " number of faces: " << nErrorFaces << endl;
1283  }
1284 
1285  return true;
1286  }
1287  else
1288  {
1289  if (debug || report)
1290  {
1291  Info<< " Face vertices OK." << endl;
1292  }
1293 
1294  return false;
1295  }
1296 }
1297 
1298 
1301  const bool report,
1302  labelHashSet* setPtr
1303 ) const
1304 {
1305  if (debug)
1306  {
1307  InfoInFunction << "Checking points" << endl;
1308  }
1309 
1310  label nFaceErrors = 0;
1311  label nCellErrors = 0;
1312 
1313  const labelListList& pf = pointFaces();
1314 
1315  forAll(pf, pointi)
1316  {
1317  if (pf[pointi].empty())
1318  {
1319  if (setPtr)
1320  {
1321  setPtr->insert(pointi);
1322  }
1323 
1324  nFaceErrors++;
1325  }
1326  }
1327 
1328 
1329  forAll(pf, pointi)
1330  {
1331  const labelList& pc = pointCells(pointi);
1332 
1333  if (pc.empty())
1334  {
1335  if (setPtr)
1336  {
1337  setPtr->insert(pointi);
1338  }
1339 
1340  nCellErrors++;
1341  }
1342  }
1343 
1344  reduce(nFaceErrors, sumOp<label>());
1345  reduce(nCellErrors, sumOp<label>());
1346 
1347  if (nFaceErrors > 0 || nCellErrors > 0)
1348  {
1349  if (debug || report)
1350  {
1351  Info<< " ***Unused points found in the mesh, "
1352  "number unused by faces: " << nFaceErrors
1353  << " number unused by cells: " << nCellErrors
1354  << endl;
1355  }
1356 
1357  return true;
1358  }
1359  else
1360  {
1361  if (debug || report)
1362  {
1363  Info<< " Point usage OK." << endl;
1364  }
1365 
1366  return false;
1367  }
1368 }
1369 
1370 
1373  const label facei,
1374  const Map<label>& nCommonPoints,
1375  label& nBaffleFaces,
1376  labelHashSet* setPtr
1377 ) const
1378 {
1379  bool error = false;
1380 
1381  forAllConstIter(Map<label>, nCommonPoints, iter)
1382  {
1383  label nbFacei = iter.key();
1384  label nCommon = iter();
1385 
1386  const face& curFace = faces()[facei];
1387  const face& nbFace = faces()[nbFacei];
1388 
1389  if (nCommon == nbFace.size() || nCommon == curFace.size())
1390  {
1391  if (nbFace.size() != curFace.size())
1392  {
1393  error = true;
1394  }
1395  else
1396  {
1397  nBaffleFaces++;
1398  }
1399 
1400  if (setPtr)
1401  {
1402  setPtr->insert(facei);
1403  setPtr->insert(nbFacei);
1404  }
1405  }
1406  }
1407 
1408  return error;
1409 }
1410 
1411 
1414  const label facei,
1415  const Map<label>& nCommonPoints,
1416  labelHashSet* setPtr
1417 ) const
1418 {
1419  bool error = false;
1420 
1421  forAllConstIter(Map<label>, nCommonPoints, iter)
1422  {
1423  label nbFacei = iter.key();
1424  label nCommon = iter();
1425 
1426  const face& curFace = faces()[facei];
1427  const face& nbFace = faces()[nbFacei];
1428 
1429  if
1430  (
1431  nCommon >= 2
1432  && nCommon != nbFace.size()
1433  && nCommon != curFace.size()
1434  )
1435  {
1436  forAll(curFace, fp)
1437  {
1438  // Get the index in the neighbouring face shared with curFace
1439  label nb = findIndex(nbFace, curFace[fp]);
1440 
1441  if (nb != -1)
1442  {
1443 
1444  // Check the whole face from nb onwards for shared vertices
1445  // with neighbouring face. Rule is that any shared vertices
1446  // should be consecutive on both faces i.e. if they are
1447  // vertices fp,fp+1,fp+2 on one face they should be
1448  // vertices nb, nb+1, nb+2 (or nb+2, nb+1, nb) on the
1449  // other face.
1450 
1451 
1452  // Vertices before and after on curFace
1453  label fpPlus1 = curFace.fcIndex(fp);
1454  label fpMin1 = curFace.rcIndex(fp);
1455 
1456  // Vertices before and after on nbFace
1457  label nbPlus1 = nbFace.fcIndex(nb);
1458  label nbMin1 = nbFace.rcIndex(nb);
1459 
1460  // Find order of walking by comparing next points on both
1461  // faces.
1462  label curInc = labelMax;
1463  label nbInc = labelMax;
1464 
1465  if (nbFace[nbPlus1] == curFace[fpPlus1])
1466  {
1467  curInc = 1;
1468  nbInc = 1;
1469  }
1470  else if (nbFace[nbPlus1] == curFace[fpMin1])
1471  {
1472  curInc = -1;
1473  nbInc = 1;
1474  }
1475  else if (nbFace[nbMin1] == curFace[fpMin1])
1476  {
1477  curInc = -1;
1478  nbInc = -1;
1479  }
1480  else
1481  {
1482  curInc = 1;
1483  nbInc = -1;
1484  }
1485 
1486 
1487  // Pass1: loop until start of common vertices found.
1488  label curNb = nb;
1489  label curFp = fp;
1490 
1491  do
1492  {
1493  curFp += curInc;
1494 
1495  if (curFp >= curFace.size())
1496  {
1497  curFp = 0;
1498  }
1499  else if (curFp < 0)
1500  {
1501  curFp = curFace.size()-1;
1502  }
1503 
1504  curNb += nbInc;
1505 
1506  if (curNb >= nbFace.size())
1507  {
1508  curNb = 0;
1509  }
1510  else if (curNb < 0)
1511  {
1512  curNb = nbFace.size()-1;
1513  }
1514  } while (curFace[curFp] == nbFace[curNb]);
1515 
1516 
1517  // Pass2: check equality walking from curFp, curNb
1518  // in opposite order.
1519 
1520  curInc = -curInc;
1521  nbInc = -nbInc;
1522 
1523  for (label commonI = 0; commonI < nCommon; commonI++)
1524  {
1525  curFp += curInc;
1526 
1527  if (curFp >= curFace.size())
1528  {
1529  curFp = 0;
1530  }
1531  else if (curFp < 0)
1532  {
1533  curFp = curFace.size()-1;
1534  }
1535 
1536  curNb += nbInc;
1537 
1538  if (curNb >= nbFace.size())
1539  {
1540  curNb = 0;
1541  }
1542  else if (curNb < 0)
1543  {
1544  curNb = nbFace.size()-1;
1545  }
1546 
1547  if (curFace[curFp] != nbFace[curNb])
1548  {
1549  if (setPtr)
1550  {
1551  setPtr->insert(facei);
1552  setPtr->insert(nbFacei);
1553  }
1554 
1555  error = true;
1556 
1557  break;
1558  }
1559  }
1560 
1561 
1562  // Done the curFace - nbFace combination.
1563  break;
1564  }
1565  }
1566  }
1567  }
1568 
1569  return error;
1570 }
1571 
1572 
1575  const bool report,
1576  labelHashSet* setPtr
1577 ) const
1578 {
1579  if (debug)
1580  {
1581  InfoInFunction << "Checking face-face connectivity" << endl;
1582  }
1583 
1584  const labelListList& pf = pointFaces();
1585 
1586  label nBaffleFaces = 0;
1587  label nErrorDuplicate = 0;
1588  label nErrorOrder = 0;
1589  Map<label> nCommonPoints(100);
1590 
1591  for (label facei = 0; facei < nFaces(); facei++)
1592  {
1593  const face& curFace = faces()[facei];
1594 
1595  // Calculate number of common points between current facei and
1596  // neighbouring face. Store on map.
1597  nCommonPoints.clear();
1598 
1599  forAll(curFace, fp)
1600  {
1601  label pointi = curFace[fp];
1602 
1603  const labelList& nbs = pf[pointi];
1604 
1605  forAll(nbs, nbI)
1606  {
1607  label nbFacei = nbs[nbI];
1608 
1609  if (facei < nbFacei)
1610  {
1611  // Only check once for each combination of two faces.
1612 
1613  Map<label>::iterator fnd = nCommonPoints.find(nbFacei);
1614 
1615  if (fnd == nCommonPoints.end())
1616  {
1617  // First common vertex found.
1618  nCommonPoints.insert(nbFacei, 1);
1619  }
1620  else
1621  {
1622  fnd()++;
1623  }
1624  }
1625  }
1626  }
1627 
1628  // Perform various checks on common points
1629 
1630  // Check all vertices shared (duplicate point)
1631  if (checkDuplicateFaces(facei, nCommonPoints, nBaffleFaces, setPtr))
1632  {
1633  nErrorDuplicate++;
1634  }
1635 
1636  // Check common vertices are consecutive on both faces
1637  if (checkCommonOrder(facei, nCommonPoints, setPtr))
1638  {
1639  nErrorOrder++;
1640  }
1641  }
1642 
1643  reduce(nBaffleFaces, sumOp<label>());
1644  reduce(nErrorDuplicate, sumOp<label>());
1645  reduce(nErrorOrder, sumOp<label>());
1646 
1647  if (nBaffleFaces)
1648  {
1649  Info<< " Number of identical duplicate faces (baffle faces): "
1650  << nBaffleFaces << endl;
1651  }
1652 
1653  if (nErrorDuplicate > 0 || nErrorOrder > 0)
1654  {
1655  // These are actually warnings, not errors.
1656  if (nErrorDuplicate > 0)
1657  {
1658  Info<< " <<Number of duplicate (not baffle) faces found: "
1659  << nErrorDuplicate
1660  << ". This might indicate a problem." << endl;
1661  }
1662 
1663  if (nErrorOrder > 0)
1664  {
1665  Info<< " <<Number of faces with non-consecutive shared points: "
1666  << nErrorOrder << ". This might indicate a problem." << endl;
1667  }
1668 
1669  return false; // return true;
1670  }
1671  else
1672  {
1673  if (debug || report)
1674  {
1675  Info<< " Face-face connectivity OK." << endl;
1676  }
1677 
1678  return false;
1679  }
1680 }
1681 
1682 
1683 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1684 
1685 bool Foam::primitiveMesh::checkClosedBoundary(const bool report) const
1686 {
1687  return checkClosedBoundary(faceAreas(), report, PackedBoolList(0));
1688 }
1689 
1690 
1693  const bool report,
1694  labelHashSet* setPtr,
1695  labelHashSet* aspectSetPtr,
1696  const Vector<label>& solutionD
1697 ) const
1698 {
1699  return checkClosedCells
1700  (
1701  faceAreas(),
1702  cellVolumes(),
1703  report,
1704  setPtr,
1705  aspectSetPtr,
1706  solutionD
1707  );
1708 }
1709 
1710 
1713  const bool report,
1714  labelHashSet* setPtr
1715 ) const
1716 {
1717  return checkFaceAreas
1718  (
1719  faceAreas(),
1720  report,
1721  false, // detailedReport,
1722  setPtr
1723  );
1724 }
1725 
1726 
1729  const bool report,
1730  labelHashSet* setPtr
1731 ) const
1732 {
1733  return checkCellVolumes
1734  (
1735  cellVolumes(),
1736  report,
1737  false, // detailedReport,
1738  setPtr
1739  );
1740 }
1741 
1742 
1745  const bool report,
1746  labelHashSet* setPtr
1747 ) const
1748 {
1749  return checkFaceOrthogonality
1750  (
1751  faceAreas(),
1752  cellCentres(),
1753  report,
1754  setPtr
1755  );
1756 }
1757 
1758 
1761  const bool report,
1762  const scalar minPyrVol,
1763  labelHashSet* setPtr
1764 ) const
1765 {
1766  return checkFacePyramids
1767  (
1768  points(),
1769  cellCentres(),
1770  report,
1771  false, // detailedReport,
1772  minPyrVol,
1773  setPtr
1774  );
1775 }
1776 
1777 
1780  const bool report,
1781  labelHashSet* setPtr
1782 ) const
1783 {
1784  return checkFaceSkewness
1785  (
1786  points(),
1787  faceCentres(),
1788  faceAreas(),
1789  cellCentres(),
1790  report,
1791  setPtr
1792  );
1793 }
1794 
1795 
1798  const bool report,
1799  const scalar maxDeg,
1800  labelHashSet* setPtr
1801 ) const
1802 {
1803  return checkFaceAngles
1804  (
1805  points(),
1806  faceAreas(),
1807  report,
1808  maxDeg,
1809  setPtr
1810  );
1811 }
1812 
1813 
1816  const bool report,
1817  const scalar warnFlatness,
1818  labelHashSet* setPtr
1819 ) const
1820 {
1821  return checkFaceFlatness
1822  (
1823  points(),
1824  faceCentres(),
1825  faceAreas(),
1826  report,
1827  warnFlatness,
1828  setPtr
1829  );
1830 }
1831 
1832 
1835  const bool report,
1836  labelHashSet* setPtr
1837 ) const
1838 {
1839  return checkConcaveCells
1840  (
1841  faceAreas(),
1842  faceCentres(),
1843  report,
1844  setPtr
1845  );
1846 }
1847 
1848 
1849 bool Foam::primitiveMesh::checkTopology(const bool report) const
1850 {
1851  label noFailedChecks = 0;
1852 
1853  if (checkPoints(report)) noFailedChecks++;
1854  if (checkUpperTriangular(report)) noFailedChecks++;
1855  if (checkCellsZipUp(report)) noFailedChecks++;
1856  if (checkFaceVertices(report)) noFailedChecks++;
1857  if (checkFaceFaces(report)) noFailedChecks++;
1858 
1859  if (noFailedChecks == 0)
1860  {
1861  if (debug || report)
1862  {
1863  Info<< " Mesh topology OK." << endl;
1864  }
1865 
1866  return false;
1867  }
1868  else
1869  {
1870  if (debug || report)
1871  {
1872  Info<< " Failed " << noFailedChecks
1873  << " mesh topology checks." << endl;
1874  }
1875 
1876  return true;
1877  }
1878 }
1879 
1880 
1881 bool Foam::primitiveMesh::checkGeometry(const bool report) const
1882 {
1883  label noFailedChecks = 0;
1884 
1885  if (checkClosedBoundary(report)) noFailedChecks++;
1886  if (checkClosedCells(report)) noFailedChecks++;
1887  if (checkFaceAreas(report)) noFailedChecks++;
1888  if (checkCellVolumes(report)) noFailedChecks++;
1889  if (checkFaceOrthogonality(report)) noFailedChecks++;
1890  if (checkFacePyramids(report)) noFailedChecks++;
1891  if (checkFaceSkewness(report)) noFailedChecks++;
1892 
1893  if (noFailedChecks == 0)
1894  {
1895  if (debug || report)
1896  {
1897  Info<< " Mesh geometry OK." << endl;
1898  }
1899  return false;
1900  }
1901  else
1902  {
1903  if (debug || report)
1904  {
1905  Info<< " Failed " << noFailedChecks
1906  << " mesh geometry checks." << endl;
1907  }
1908 
1909  return true;
1910  }
1911 }
1912 
1913 
1914 bool Foam::primitiveMesh::checkMesh(const bool report) const
1915 {
1916  if (debug)
1917  {
1918  InfoInFunction << "Checking primitiveMesh" << endl;
1919  }
1920 
1921  label noFailedChecks = checkTopology(report) + checkGeometry(report);
1922 
1923  if (noFailedChecks == 0)
1924  {
1925  if (debug || report)
1926  {
1927  Info<< "Mesh OK." << endl;
1928  }
1929 
1930  return false;
1931  }
1932  else
1933  {
1934  if (debug || report)
1935  {
1936  Info<< " Failed " << noFailedChecks
1937  << " mesh checks." << endl;
1938  }
1939 
1940  return true;
1941  }
1942 }
1943 
1944 
1945 Foam::scalar Foam::primitiveMesh::setClosedThreshold(const scalar val)
1946 {
1947  scalar prev = closedThreshold_;
1948  closedThreshold_ = val;
1949 
1950  return prev;
1951 }
1952 
1953 
1954 Foam::scalar Foam::primitiveMesh::setAspectThreshold(const scalar val)
1955 {
1956  scalar prev = aspectThreshold_;
1957  aspectThreshold_ = val;
1958 
1959  return prev;
1960 }
1961 
1962 
1963 Foam::scalar Foam::primitiveMesh::setNonOrthThreshold(const scalar val)
1964 {
1965  scalar prev = nonOrthThreshold_;
1966  nonOrthThreshold_ = val;
1967 
1968  return prev;
1969 }
1970 
1971 
1972 Foam::scalar Foam::primitiveMesh::setSkewThreshold(const scalar val)
1973 {
1974  scalar prev = skewThreshold_;
1975  skewThreshold_ = val;
1976 
1977  return prev;
1978 }
1979 
1980 
1981 // ************************************************************************* //
bool checkFaceAreas(const vectorField &faceAreas, const bool report, const bool detailedReport, labelHashSet *setPtr) const
Check for negative face areas.
static scalar aspectThreshold_
Aspect ratio warning threshold.
dimensionedScalar acos(const dimensionedScalar &ds)
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
bool checkFaceAngles(const pointField &points, const vectorField &faceAreas, const bool report, const scalar maxDeg, labelHashSet *setPtr) const
Check face angles.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:325
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:112
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
Definition: UListI.H:65
A list that is sorted upon construction or when explicitly requested with the sort() method...
Definition: List.H:80
Unit conversion functions.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
static scalar setClosedThreshold(const scalar)
Set the closedness ratio warning threshold.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual const pointField & points() const =0
Return mesh points.
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
const dimensionedScalar c
Speed of light in a vacuum.
dimensionedScalar asin(const dimensionedScalar &ds)
static scalar setAspectThreshold(const scalar)
Set the aspect ratio warning threshold.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Various functions to operate on Lists.
static scalar setSkewThreshold(const scalar)
Set the skewness warning threshold as percentage.
Class to handle errors and exceptions in a simple, consistent stream-based manner.
Definition: error.H:67
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
bool checkFaceSkewness(const pointField &points, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs, const bool report, labelHashSet *setPtr) const
Check face skewness.
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
bool checkClosedCells(const vectorField &faceAreas, const scalarField &cellVolumes, const bool report, labelHashSet *setPtr, labelHashSet *aspectSetPtr, const Vector< label > &meshD) const
Check cells for closedness.
bool checkFacePyramids(const pointField &points, const vectorField &ctrs, const bool report, const bool detailedReport, const scalar minPyrVol, labelHashSet *setPtr) const
Check face pyramid volume.
Type gSum(const FieldField< Field, Type > &f)
dimensionedScalar cos(const dimensionedScalar &ds)
const cellShapeList & cells
bool checkFaceOrthogonality(const vectorField &fAreas, const vectorField &cellCtrs, const bool report, labelHashSet *setPtr) const
Check for non-orthogonality.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
static scalar closedThreshold_
Static data to control mesh checking.
label nPoints
void clear()
Clear all entries from table.
Definition: HashTable.C:468
dimensionSet cmptMag(const dimensionSet &)
Definition: dimensionSet.C:275
virtual bool checkMesh(const bool report=false) const
Check mesh for correctness. Returns false for no error.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
static const label labelMax
Definition: label.H:62
static const zero Zero
Definition: zero.H:97
static scalar setNonOrthThreshold(const scalar)
Set the non-orthogonality warning threshold in degrees.
const vectorField & cellCentres() const
bool checkCellVolumes(const scalarField &vols, const bool report, const bool detailedReport, labelHashSet *setPtr) const
Check for negative cell volumes.
static scalar planarCosAngle_
Threshold where faces are considered coplanar.
dimensionedScalar sin(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:260
labelList f(nPoints)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
virtual bool checkFaceVertices(const bool report=false, labelHashSet *setPtr=nullptr) const
Check uniqueness of face vertices.
const vectorField & faceCentres() const
A bit-packed bool list.
virtual bool checkTopology(const bool report=false) const
Check mesh topology for correctness.
bool checkConcaveCells(const vectorField &fAreas, const pointField &fCentres, const bool report, labelHashSet *setPtr) const
Check for concave cells by the planes of faces.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
bool checkCommonOrder(const label, const Map< label > &, labelHashSet *) const
Check that shared points are in consecutive order.
const vectorField & faceAreas() const
bool checkDuplicateFaces(const label, const Map< label > &, label &nBaffleFaces, labelHashSet *) const
Check if all points on face are shared with another face.
messageStream Info
virtual bool checkCellsZipUp(const bool report=false, labelHashSet *setPtr=nullptr) const
Check cell zip-up.
virtual bool checkUpperTriangular(const bool report=false, labelHashSet *setPtr=nullptr) const
Check face ordering.
virtual bool checkPoints(const bool report=false, labelHashSet *setPtr=nullptr) const
Check for unused points.
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual bool checkGeometry(const bool report=false) const
Check mesh geometry (& implicitly topology) for correctness.
bool checkFaceFlatness(const pointField &points, const vectorField &faceCentres, const vectorField &faceAreas, const bool report, const scalar warnFlatness, labelHashSet *setPtr) const
Check face warpage.
label size() const
Number of entries.
Definition: PackedListI.H:711
virtual bool checkFaceFaces(const bool report=false, labelHashSet *setPtr=nullptr) const
Check face-face connectivity.
static scalar skewThreshold_
Skewness warning threshold.
static scalar nonOrthThreshold_
Non-orthogonality warning threshold in deg.
A class for managing temporary objects.
Definition: PtrList.H:53
bool checkClosedBoundary(const vectorField &, const bool, const PackedBoolList &) const
Check boundary for closedness.
const scalarField & cellVolumes() const
#define InfoInFunction
Report an information message using Foam::Info.