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