polyMeshCheck.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-2021 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 "polyMeshCheck.H"
28 #include "pyramidPointFaceRef.H"
29 #include "tetPointRef.H"
30 #include "syncTools.H"
31 #include "unitConversion.H"
32 #include "primitiveMeshTools.H"
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 namespace polyMeshCheck
39 {
40 
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 
43 scalar checkNonOrtho
44 (
45  const polyMesh& mesh,
46  const bool report,
47  const scalar severeNonorthogonalityThreshold,
48  const label facei,
49  const vector& s, // face area vector
50  const vector& d, // cc-cc vector
51 
52  label& severeNonOrth,
53  label& errorNonOrth,
54  labelHashSet* setPtr
55 )
56 {
57  scalar dDotS = (d & s)/(mag(d)*mag(s) + vSmall);
58 
59  if (dDotS < severeNonorthogonalityThreshold)
60  {
61  label nei = -1;
62 
63  if (mesh.isInternalFace(facei))
64  {
65  nei = mesh.faceNeighbour()[facei];
66  }
67 
68  if (dDotS > small)
69  {
70  if (report)
71  {
72  // Severe non-orthogonality but mesh still OK
73  Pout<< "Severe non-orthogonality for face " << facei
74  << " between cells " << mesh.faceOwner()[facei]
75  << " and " << nei
76  << ": Angle = "
77  << radToDeg(::acos(dDotS))
78  << " deg." << endl;
79  }
80 
81  severeNonOrth++;
82  }
83  else
84  {
85  // Non-orthogonality greater than 90 deg
86  if (report)
87  {
89  << "Severe non-orthogonality detected for face "
90  << facei
91  << " between cells " << mesh.faceOwner()[facei]
92  << " and " << nei
93  << ": Angle = "
94  << radToDeg(::acos(dDotS))
95  << " deg." << endl;
96  }
97 
98  errorNonOrth++;
99  }
100 
101  if (setPtr)
102  {
103  setPtr->insert(facei);
104  }
105  }
106  return dDotS;
107 }
108 
109 
110 bool checkFaceTet
111 (
112  const polyMesh& mesh,
113  const bool report,
114  const scalar minTetQuality,
115  const pointField& p,
116  const label facei,
117  const point& fc, // face centre
118  const point& cc, // cell centre
119 
120  labelHashSet* setPtr
121 )
122 {
123  const face& f = mesh.faces()[facei];
124 
125  forAll(f, fp)
126  {
127  scalar tetQual = tetPointRef
128  (
129  p[f[fp]],
130  p[f.nextLabel(fp)],
131  fc,
132  cc
133  ).quality();
134 
135  if (tetQual < minTetQuality)
136  {
137  if (report)
138  {
139  Pout<< "bool polyMeshCheck::checkFaceTets("
140  << "const bool, const scalar, const pointField&"
141  << ", const pointField&"
142  << ", const labelList&, labelHashSet*) : "
143  << "face " << facei
144  << " has a triangle that points the wrong way."
145  << endl
146  << "Tet quality: " << tetQual
147  << " Face " << facei
148  << endl;
149  }
150 
151  if (setPtr)
152  {
153  setPtr->insert(facei);
154  }
155  return true;
156  }
157  }
158  return false;
159 }
160 
161 
163 (
164  const polyMesh& mesh,
165  const labelList& changedFaces
166 )
167 {
168  const labelList& own = mesh.faceOwner();
169  const labelList& nei = mesh.faceNeighbour();
170 
171  labelHashSet affectedCells(2*changedFaces.size());
172 
173  forAll(changedFaces, i)
174  {
175  label facei = changedFaces[i];
176 
177  affectedCells.insert(own[facei]);
178 
179  if (mesh.isInternalFace(facei))
180  {
181  affectedCells.insert(nei[facei]);
182  }
183  }
184  return affectedCells.toc();
185 }
186 
187 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
188 
189 } // End namespace polyMeshCheck
190 } // End namespace Foam
191 
192 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
193 
195 (
196  const bool report,
197  const scalar orthWarn,
198  const polyMesh& mesh,
199  const vectorField& cellCentres,
200  const vectorField& faceAreas,
201  const labelList& checkFaces,
202  const List<labelPair>& baffles,
203  labelHashSet* setPtr
204 )
205 {
206  // for all internal and coupled faces check theat the d dot S product
207  // is positive
208 
209  const labelList& own = mesh.faceOwner();
210  const labelList& nei = mesh.faceNeighbour();
211  const polyBoundaryMesh& patches = mesh.boundaryMesh();
212 
213  // Severe nonorthogonality threshold
214  const scalar severeNonorthogonalityThreshold = ::cos(degToRad(orthWarn));
215 
216  // Calculate coupled cell centre
217  pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
218 
219  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
220  {
221  neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
222  }
223 
225 
226  scalar minDDotS = great;
227 
228  scalar sumDDotS = 0;
229  label nDDotS = 0;
230 
231  label severeNonOrth = 0;
232 
233  label errorNonOrth = 0;
234 
235  forAll(checkFaces, i)
236  {
237  label facei = checkFaces[i];
238 
239  const point& ownCc = cellCentres[own[facei]];
240 
241  if (mesh.isInternalFace(facei))
242  {
243  scalar dDotS = checkNonOrtho
244  (
245  mesh,
246  report,
247  severeNonorthogonalityThreshold,
248  facei,
249  faceAreas[facei],
250  cellCentres[nei[facei]] - ownCc,
251 
252  severeNonOrth,
253  errorNonOrth,
254  setPtr
255  );
256 
257  if (dDotS < minDDotS)
258  {
259  minDDotS = dDotS;
260  }
261 
262  sumDDotS += dDotS;
263  nDDotS++;
264  }
265  else
266  {
267  label patchi = patches.whichPatch(facei);
268 
269  if (patches[patchi].coupled())
270  {
271  scalar dDotS = checkNonOrtho
272  (
273  mesh,
274  report,
275  severeNonorthogonalityThreshold,
276  facei,
277  faceAreas[facei],
278  neiCc[facei-mesh.nInternalFaces()] - ownCc,
279 
280  severeNonOrth,
281  errorNonOrth,
282  setPtr
283  );
284 
285  if (dDotS < minDDotS)
286  {
287  minDDotS = dDotS;
288  }
289 
290  sumDDotS += dDotS;
291  nDDotS++;
292  }
293  }
294  }
295 
296  forAll(baffles, i)
297  {
298  label face0 = baffles[i].first();
299  label face1 = baffles[i].second();
300 
301  const point& ownCc = cellCentres[own[face0]];
302 
303  scalar dDotS = checkNonOrtho
304  (
305  mesh,
306  report,
307  severeNonorthogonalityThreshold,
308  face0,
309  faceAreas[face0],
310  cellCentres[own[face1]] - ownCc,
311 
312  severeNonOrth,
313  errorNonOrth,
314  setPtr
315  );
316 
317  if (dDotS < minDDotS)
318  {
319  minDDotS = dDotS;
320  }
321 
322  sumDDotS += dDotS;
323  nDDotS++;
324  }
325 
326  reduce(minDDotS, minOp<scalar>());
327  reduce(sumDDotS, sumOp<scalar>());
328  reduce(nDDotS, sumOp<label>());
329  reduce(severeNonOrth, sumOp<label>());
330  reduce(errorNonOrth, sumOp<label>());
331 
332  // Only report if there are some internal faces
333  if (nDDotS > 0)
334  {
335  if (report && minDDotS < severeNonorthogonalityThreshold)
336  {
337  Info<< "Number of non-orthogonality errors: " << errorNonOrth
338  << ". Number of severely non-orthogonal faces: "
339  << severeNonOrth << "." << endl;
340  }
341  }
342 
343  if (report)
344  {
345  if (nDDotS > 0)
346  {
347  Info<< "Mesh non-orthogonality Max: "
348  << radToDeg(::acos(minDDotS))
349  << " average: " << radToDeg(::acos(sumDDotS/nDDotS))
350  << endl;
351  }
352  }
353 
354  if (errorNonOrth > 0)
355  {
356  if (report)
357  {
359  << "Error in non-orthogonality detected" << endl;
360  }
361 
362  return true;
363  }
364  else
365  {
366  if (report)
367  {
368  Info<< "Non-orthogonality check OK.\n" << endl;
369  }
370 
371  return false;
372  }
373 }
374 
375 
377 (
378  const bool report,
379  const scalar minPyrVol,
380  const polyMesh& mesh,
381  const vectorField& cellCentres,
382  const pointField& p,
383  const labelList& checkFaces,
384  const List<labelPair>& baffles,
385  labelHashSet* setPtr
386 )
387 {
388  // check whether face area vector points to the cell with higher label
389  const labelList& own = mesh.faceOwner();
390  const labelList& nei = mesh.faceNeighbour();
391 
392  const faceList& f = mesh.faces();
393 
394  label nErrorPyrs = 0;
395 
396  forAll(checkFaces, i)
397  {
398  label facei = checkFaces[i];
399 
400  // Create the owner pyramid - it will have negative volume
401  scalar pyrVol = pyramidPointFaceRef
402  (
403  f[facei],
404  cellCentres[own[facei]]
405  ).mag(p);
406 
407  if (pyrVol > -minPyrVol)
408  {
409  if (report)
410  {
411  Pout<< "bool polyMeshCheck::checkFacePyramids("
412  << "const bool, const scalar, const pointField&"
413  << ", const labelList&, labelHashSet*): "
414  << "face " << facei << " points the wrong way. " << endl
415  << "Pyramid volume: " << -pyrVol
416  << " Face " << f[facei] << " area: " << f[facei].mag(p)
417  << " Owner cell: " << own[facei] << endl
418  << "Owner cell vertex labels: "
419  << mesh.cells()[own[facei]].labels(f)
420  << endl;
421  }
422 
423 
424  if (setPtr)
425  {
426  setPtr->insert(facei);
427  }
428 
429  nErrorPyrs++;
430  }
431 
432  if (mesh.isInternalFace(facei))
433  {
434  // Create the neighbour pyramid - it will have positive volume
435  scalar pyrVol =
436  pyramidPointFaceRef(f[facei], cellCentres[nei[facei]]).mag(p);
437 
438  if (pyrVol < minPyrVol)
439  {
440  if (report)
441  {
442  Pout<< "bool polyMeshCheck::checkFacePyramids("
443  << "const bool, const scalar, const pointField&"
444  << ", const labelList&, labelHashSet*): "
445  << "face " << facei << " points the wrong way. " << endl
446  << "Pyramid volume: " << -pyrVol
447  << " Face " << f[facei] << " area: " << f[facei].mag(p)
448  << " Neighbour cell: " << nei[facei] << endl
449  << "Neighbour cell vertex labels: "
450  << mesh.cells()[nei[facei]].labels(f)
451  << endl;
452  }
453 
454  if (setPtr)
455  {
456  setPtr->insert(facei);
457  }
458 
459  nErrorPyrs++;
460  }
461  }
462  }
463 
464  forAll(baffles, i)
465  {
466  label face0 = baffles[i].first();
467  label face1 = baffles[i].second();
468 
469  const point& ownCc = cellCentres[own[face0]];
470 
471  // Create the owner pyramid - it will have negative volume
472  scalar pyrVolOwn = pyramidPointFaceRef
473  (
474  f[face0],
475  ownCc
476  ).mag(p);
477 
478  if (pyrVolOwn > -minPyrVol)
479  {
480  if (report)
481  {
482  Pout<< "bool polyMeshCheck::checkFacePyramids("
483  << "const bool, const scalar, const pointField&"
484  << ", const labelList&, labelHashSet*): "
485  << "face " << face0 << " points the wrong way. " << endl
486  << "Pyramid volume: " << -pyrVolOwn
487  << " Face " << f[face0] << " area: " << f[face0].mag(p)
488  << " Owner cell: " << own[face0] << endl
489  << "Owner cell vertex labels: "
490  << mesh.cells()[own[face0]].labels(f)
491  << endl;
492  }
493 
494 
495  if (setPtr)
496  {
497  setPtr->insert(face0);
498  }
499 
500  nErrorPyrs++;
501  }
502 
503  // Create the neighbour pyramid - it will have positive volume
504  scalar pyrVolNbr =
505  pyramidPointFaceRef(f[face0], cellCentres[own[face1]]).mag(p);
506 
507  if (pyrVolNbr < minPyrVol)
508  {
509  if (report)
510  {
511  Pout<< "bool polyMeshCheck::checkFacePyramids("
512  << "const bool, const scalar, const pointField&"
513  << ", const labelList&, labelHashSet*): "
514  << "face " << face0 << " points the wrong way. " << endl
515  << "Pyramid volume: " << -pyrVolNbr
516  << " Face " << f[face0] << " area: " << f[face0].mag(p)
517  << " Neighbour cell: " << own[face1] << endl
518  << "Neighbour cell vertex labels: "
519  << mesh.cells()[own[face1]].labels(f)
520  << endl;
521  }
522 
523  if (setPtr)
524  {
525  setPtr->insert(face0);
526  }
527 
528  nErrorPyrs++;
529  }
530  }
531 
532  reduce(nErrorPyrs, sumOp<label>());
533 
534  if (nErrorPyrs > 0)
535  {
536  if (report)
537  {
539  << "Error in face pyramids: faces pointing the wrong way."
540  << endl;
541  }
542 
543  return true;
544  }
545  else
546  {
547  if (report)
548  {
549  Info<< "Face pyramids OK.\n" << endl;
550  }
551 
552  return false;
553  }
554 }
555 
556 
558 (
559  const bool report,
560  const scalar minTetQuality,
561  const polyMesh& mesh,
562  const vectorField& cellCentres,
563  const vectorField& faceCentres,
564  const pointField& p,
565  const labelList& checkFaces,
566  const List<labelPair>& baffles,
567  labelHashSet* setPtr
568 )
569 {
570  // check whether decomposing each cell into tets results in
571  // positive volume, non-flat tets
572  const labelList& own = mesh.faceOwner();
573  const labelList& nei = mesh.faceNeighbour();
574  const polyBoundaryMesh& patches = mesh.boundaryMesh();
575 
576  // Calculate coupled cell centre
577  pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
578 
579  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
580  {
581  neiCc[facei - mesh.nInternalFaces()] = cellCentres[own[facei]];
582  }
583 
585 
586  label nErrorTets = 0;
587 
588  forAll(checkFaces, i)
589  {
590  label facei = checkFaces[i];
591 
592  // Create the owner pyramid - note: exchange cell and face centre
593  // to get positive volume.
594  bool tetError = checkFaceTet
595  (
596  mesh,
597  report,
598  minTetQuality,
599  p,
600  facei,
601  cellCentres[own[facei]], // face centre
602  faceCentres[facei], // cell centre
603  setPtr
604  );
605 
606  if (tetError)
607  {
608  nErrorTets++;
609  }
610 
611  if (mesh.isInternalFace(facei))
612  {
613  // Create the neighbour tets - they will have positive volume
614  bool tetError = checkFaceTet
615  (
616  mesh,
617  report,
618  minTetQuality,
619  p,
620  facei,
621  faceCentres[facei], // face centre
622  cellCentres[nei[facei]], // cell centre
623  setPtr
624  );
625 
626  if (tetError)
627  {
628  nErrorTets++;
629  }
630 
631  if
632  (
634  (
635  mesh,
636  facei,
637  minTetQuality,
638  report
639  ) == -1
640  )
641  {
642  if (setPtr)
643  {
644  setPtr->insert(facei);
645  }
646 
647  nErrorTets++;
648  }
649  }
650  else
651  {
652  label patchi = patches.whichPatch(facei);
653 
654  if (patches[patchi].coupled())
655  {
656  if
657  (
659  (
660  mesh,
661  facei,
662  neiCc[facei - mesh.nInternalFaces()],
663  minTetQuality,
664  report
665  ) == -1
666  )
667  {
668  if (setPtr)
669  {
670  setPtr->insert(facei);
671  }
672 
673  nErrorTets++;
674  }
675  }
676  else
677  {
678  if
679  (
681  (
682  mesh,
683  facei,
684  minTetQuality,
685  report
686  ) == -1
687  )
688  {
689  if (setPtr)
690  {
691  setPtr->insert(facei);
692  }
693 
694  nErrorTets++;
695  }
696  }
697  }
698  }
699 
700  forAll(baffles, i)
701  {
702  label face0 = baffles[i].first();
703  label face1 = baffles[i].second();
704 
705  bool tetError = checkFaceTet
706  (
707  mesh,
708  report,
709  minTetQuality,
710  p,
711  face0,
712  cellCentres[own[face0]], // face centre
713  faceCentres[face0], // cell centre
714  setPtr
715  );
716 
717  if (tetError)
718  {
719  nErrorTets++;
720  }
721 
722  // Create the neighbour tets - they will have positive volume
723  tetError = checkFaceTet
724  (
725  mesh,
726  report,
727  minTetQuality,
728  p,
729  face0,
730  faceCentres[face0], // face centre
731  cellCentres[own[face1]], // cell centre
732  setPtr
733  );
734 
735  if (tetError)
736  {
737  nErrorTets++;
738  }
739 
740  if
741  (
743  (
744  mesh,
745  face0,
746  cellCentres[own[face1]],
747  minTetQuality,
748  report
749  ) == -1
750  )
751  {
752  if (setPtr)
753  {
754  setPtr->insert(face0);
755  }
756 
757  nErrorTets++;
758  }
759  }
760 
761  reduce(nErrorTets, sumOp<label>());
762 
763  if (nErrorTets > 0)
764  {
765  if (report)
766  {
768  << "Error in face decomposition: negative tets."
769  << endl;
770  }
771 
772  return true;
773  }
774  else
775  {
776  if (report)
777  {
778  Info<< "Face tets OK.\n" << endl;
779  }
780 
781  return false;
782  }
783 }
784 
785 
787 (
788  const bool report,
789  const scalar internalSkew,
790  const scalar boundarySkew,
791  const polyMesh& mesh,
792  const pointField& points,
793  const vectorField& cellCentres,
794  const vectorField& faceCentres,
795  const vectorField& faceAreas,
796  const labelList& checkFaces,
797  const List<labelPair>& baffles,
798  labelHashSet* setPtr
799 )
800 {
801  // Warn if the skew correction vector is more than skew times
802  // larger than the face area vector
803 
804  const labelList& own = mesh.faceOwner();
805  const labelList& nei = mesh.faceNeighbour();
806  const polyBoundaryMesh& patches = mesh.boundaryMesh();
807 
808  // Calculate coupled cell centre
809  pointField neiCc;
810  syncTools::swapBoundaryCellPositions(mesh, cellCentres, neiCc);
811 
812  scalar maxSkew = 0;
813 
814  label nWarnSkew = 0;
815 
816  forAll(checkFaces, i)
817  {
818  label facei = checkFaces[i];
819 
820  if (mesh.isInternalFace(facei))
821  {
822  scalar skewness = primitiveMeshTools::faceSkewness
823  (
824  mesh,
825  points,
826  faceCentres,
827  faceAreas,
828 
829  facei,
830  cellCentres[own[facei]],
831  cellCentres[nei[facei]]
832  );
833 
834  // Check if the skewness vector is greater than the PN vector.
835  // This does not cause trouble but is a good indication of a poor
836  // mesh.
837  if (skewness > internalSkew)
838  {
839  if (report)
840  {
841  Pout<< "Severe skewness for face " << facei
842  << " skewness = " << skewness << endl;
843  }
844 
845  if (setPtr)
846  {
847  setPtr->insert(facei);
848  }
849 
850  nWarnSkew++;
851  }
852 
853  maxSkew = max(maxSkew, skewness);
854  }
855  else if (patches[patches.whichPatch(facei)].coupled())
856  {
857  scalar skewness = primitiveMeshTools::faceSkewness
858  (
859  mesh,
860  points,
861  faceCentres,
862  faceAreas,
863 
864  facei,
865  cellCentres[own[facei]],
866  neiCc[facei-mesh.nInternalFaces()]
867  );
868 
869  // Check if the skewness vector is greater than the PN vector.
870  // This does not cause trouble but is a good indication of a poor
871  // mesh.
872  if (skewness > internalSkew)
873  {
874  if (report)
875  {
876  Pout<< "Severe skewness for coupled face " << facei
877  << " skewness = " << skewness << endl;
878  }
879 
880  if (setPtr)
881  {
882  setPtr->insert(facei);
883  }
884 
885  nWarnSkew++;
886  }
887 
888  maxSkew = max(maxSkew, skewness);
889  }
890  else
891  {
893  (
894  mesh,
895  points,
896  faceCentres,
897  faceAreas,
898 
899  facei,
900  cellCentres[own[facei]]
901  );
902 
903 
904  // Check if the skewness vector is greater than the PN vector.
905  // This does not cause trouble but is a good indication of a poor
906  // mesh.
907  if (skewness > boundarySkew)
908  {
909  if (report)
910  {
911  Pout<< "Severe skewness for boundary face " << facei
912  << " skewness = " << skewness << endl;
913  }
914 
915  if (setPtr)
916  {
917  setPtr->insert(facei);
918  }
919 
920  nWarnSkew++;
921  }
922 
923  maxSkew = max(maxSkew, skewness);
924  }
925  }
926 
927  forAll(baffles, i)
928  {
929  label face0 = baffles[i].first();
930  label face1 = baffles[i].second();
931 
932  const point& ownCc = cellCentres[own[face0]];
933  const point& neiCc = cellCentres[own[face1]];
934 
935  scalar skewness = primitiveMeshTools::faceSkewness
936  (
937  mesh,
938  points,
939  faceCentres,
940  faceAreas,
941 
942  face0,
943  ownCc,
944  neiCc
945  );
946 
947  // Check if the skewness vector is greater than the PN vector.
948  // This does not cause trouble but is a good indication of a poor
949  // mesh.
950  if (skewness > internalSkew)
951  {
952  if (report)
953  {
954  Pout<< "Severe skewness for face " << face0
955  << " skewness = " << skewness << endl;
956  }
957 
958  if (setPtr)
959  {
960  setPtr->insert(face0);
961  }
962 
963  nWarnSkew++;
964  }
965 
966  maxSkew = max(maxSkew, skewness);
967  }
968 
969 
970  reduce(maxSkew, maxOp<scalar>());
971  reduce(nWarnSkew, sumOp<label>());
972 
973  if (nWarnSkew > 0)
974  {
975  if (report)
976  {
978  << 100*maxSkew
979  << " percent.\nThis may impair the quality of the result." << nl
980  << nWarnSkew << " highly skew faces detected."
981  << endl;
982  }
983 
984  return true;
985  }
986  else
987  {
988  if (report)
989  {
990  Info<< "Max skewness = " << 100*maxSkew
991  << " percent. Face skewness OK.\n" << endl;
992  }
993 
994  return false;
995  }
996 }
997 
998 
1001  const bool report,
1002  const scalar warnWeight,
1003  const polyMesh& mesh,
1004  const vectorField& cellCentres,
1005  const vectorField& faceCentres,
1006  const vectorField& faceAreas,
1007  const labelList& checkFaces,
1008  const List<labelPair>& baffles,
1009  labelHashSet* setPtr
1010 )
1011 {
1012  // Warn if the delta factor (0..1) is too large.
1013 
1014  const labelList& own = mesh.faceOwner();
1015  const labelList& nei = mesh.faceNeighbour();
1016  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1017 
1018  // Calculate coupled cell centre
1019  pointField neiCc(mesh.nFaces()-mesh.nInternalFaces());
1020 
1021  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
1022  {
1023  neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
1024  }
1026 
1027 
1028  scalar minWeight = great;
1029 
1030  label nWarnWeight = 0;
1031 
1032  forAll(checkFaces, i)
1033  {
1034  label facei = checkFaces[i];
1035 
1036  const point& fc = faceCentres[facei];
1037  const vector& fa = faceAreas[facei];
1038 
1039  scalar dOwn = mag(fa & (fc-cellCentres[own[facei]]));
1040 
1041  if (mesh.isInternalFace(facei))
1042  {
1043  scalar dNei = mag(fa & (cellCentres[nei[facei]]-fc));
1044  scalar weight = min(dNei,dOwn)/(dNei+dOwn+vSmall);
1045 
1046  if (weight < warnWeight)
1047  {
1048  if (report)
1049  {
1050  Pout<< "Small weighting factor for face " << facei
1051  << " weight = " << weight << endl;
1052  }
1053 
1054  if (setPtr)
1055  {
1056  setPtr->insert(facei);
1057  }
1058 
1059  nWarnWeight++;
1060  }
1061 
1062  minWeight = min(minWeight, weight);
1063  }
1064  else
1065  {
1066  label patchi = patches.whichPatch(facei);
1067 
1068  if (patches[patchi].coupled())
1069  {
1070  scalar dNei = mag(fa & (neiCc[facei-mesh.nInternalFaces()]-fc));
1071  scalar weight = min(dNei,dOwn)/(dNei+dOwn+vSmall);
1072 
1073  if (weight < warnWeight)
1074  {
1075  if (report)
1076  {
1077  Pout<< "Small weighting factor for face " << facei
1078  << " weight = " << weight << endl;
1079  }
1080 
1081  if (setPtr)
1082  {
1083  setPtr->insert(facei);
1084  }
1085 
1086  nWarnWeight++;
1087  }
1088 
1089  minWeight = min(minWeight, weight);
1090  }
1091  }
1092  }
1093 
1094  forAll(baffles, i)
1095  {
1096  label face0 = baffles[i].first();
1097  label face1 = baffles[i].second();
1098 
1099  const point& ownCc = cellCentres[own[face0]];
1100  const point& fc = faceCentres[face0];
1101  const vector& fa = faceAreas[face0];
1102 
1103  scalar dOwn = mag(fa & (fc-ownCc));
1104  scalar dNei = mag(fa & (cellCentres[own[face1]]-fc));
1105  scalar weight = min(dNei,dOwn)/(dNei+dOwn+vSmall);
1106 
1107  if (weight < warnWeight)
1108  {
1109  if (report)
1110  {
1111  Pout<< "Small weighting factor for face " << face0
1112  << " weight = " << weight << endl;
1113  }
1114 
1115  if (setPtr)
1116  {
1117  setPtr->insert(face0);
1118  }
1119 
1120  nWarnWeight++;
1121  }
1122 
1123  minWeight = min(minWeight, weight);
1124  }
1125 
1126  reduce(minWeight, minOp<scalar>());
1127  reduce(nWarnWeight, sumOp<label>());
1128 
1129  if (minWeight < warnWeight)
1130  {
1131  if (report)
1132  {
1134  << minWeight << '.' << nl
1135  << nWarnWeight << " faces with small weights detected."
1136  << endl;
1137  }
1138 
1139  return true;
1140  }
1141  else
1142  {
1143  if (report)
1144  {
1145  Info<< "Min weight = " << minWeight
1146  << ". Weights OK.\n" << endl;
1147  }
1148 
1149  return false;
1150  }
1151 }
1152 
1153 
1156  const bool report,
1157  const scalar warnRatio,
1158  const polyMesh& mesh,
1159  const scalarField& cellVolumes,
1160  const labelList& checkFaces,
1161  const List<labelPair>& baffles,
1162  labelHashSet* setPtr
1163 )
1164 {
1165  // Warn if the volume ratio between neighbouring cells is too large
1166 
1167  const labelList& own = mesh.faceOwner();
1168  const labelList& nei = mesh.faceNeighbour();
1169  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1170 
1171  // Calculate coupled cell vol
1172  scalarField neiVols(mesh.nFaces()-mesh.nInternalFaces());
1173 
1174  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
1175  {
1176  neiVols[facei-mesh.nInternalFaces()] = cellVolumes[own[facei]];
1177  }
1178  syncTools::swapBoundaryFaceList(mesh, neiVols);
1179 
1180 
1181  scalar minRatio = great;
1182 
1183  label nWarnRatio = 0;
1184 
1185  forAll(checkFaces, i)
1186  {
1187  label facei = checkFaces[i];
1188 
1189  scalar ownVol = mag(cellVolumes[own[facei]]);
1190 
1191  scalar neiVol = -great;
1192 
1193  if (mesh.isInternalFace(facei))
1194  {
1195  neiVol = mag(cellVolumes[nei[facei]]);
1196  }
1197  else
1198  {
1199  label patchi = patches.whichPatch(facei);
1200 
1201  if (patches[patchi].coupled())
1202  {
1203  neiVol = mag(neiVols[facei-mesh.nInternalFaces()]);
1204  }
1205  }
1206 
1207  if (neiVol >= 0)
1208  {
1209  scalar ratio = min(ownVol, neiVol) / (max(ownVol, neiVol) + vSmall);
1210 
1211  if (ratio < warnRatio)
1212  {
1213  if (report)
1214  {
1215  Pout<< "Small ratio for face " << facei
1216  << " ratio = " << ratio << endl;
1217  }
1218 
1219  if (setPtr)
1220  {
1221  setPtr->insert(facei);
1222  }
1223 
1224  nWarnRatio++;
1225  }
1226 
1227  minRatio = min(minRatio, ratio);
1228  }
1229  }
1230 
1231  forAll(baffles, i)
1232  {
1233  label face0 = baffles[i].first();
1234  label face1 = baffles[i].second();
1235 
1236  scalar ownVol = mag(cellVolumes[own[face0]]);
1237 
1238  scalar neiVol = mag(cellVolumes[own[face1]]);
1239 
1240  if (neiVol >= 0)
1241  {
1242  scalar ratio = min(ownVol, neiVol) / (max(ownVol, neiVol) + vSmall);
1243 
1244  if (ratio < warnRatio)
1245  {
1246  if (report)
1247  {
1248  Pout<< "Small ratio for face " << face0
1249  << " ratio = " << ratio << endl;
1250  }
1251 
1252  if (setPtr)
1253  {
1254  setPtr->insert(face0);
1255  }
1256 
1257  nWarnRatio++;
1258  }
1259 
1260  minRatio = min(minRatio, ratio);
1261  }
1262  }
1263 
1264  reduce(minRatio, minOp<scalar>());
1265  reduce(nWarnRatio, sumOp<label>());
1266 
1267  if (minRatio < warnRatio)
1268  {
1269  if (report)
1270  {
1272  << minRatio << '.' << nl
1273  << nWarnRatio << " faces with small ratios detected."
1274  << endl;
1275  }
1276 
1277  return true;
1278  }
1279  else
1280  {
1281  if (report)
1282  {
1283  Info<< "Min ratio = " << minRatio
1284  << ". Ratios OK.\n" << endl;
1285  }
1286 
1287  return false;
1288  }
1289 }
1290 
1291 
1294  const bool report,
1295  const scalar maxDeg,
1296  const polyMesh& mesh,
1297  const vectorField& faceAreas,
1298  const pointField& p,
1299  const labelList& checkFaces,
1300  labelHashSet* setPtr
1301 )
1302 {
1303  if (maxDeg < -small || maxDeg > 180+small)
1304  {
1306  << "maxDeg should be [0..180] but is now " << maxDeg
1307  << abort(FatalError);
1308  }
1309 
1310  const scalar maxSin = Foam::sin(degToRad(maxDeg));
1311 
1312  const faceList& fcs = mesh.faces();
1313 
1314  scalar maxEdgeSin = 0.0;
1315 
1316  label nConcave = 0;
1317 
1318  label errorFacei = -1;
1319 
1320  forAll(checkFaces, i)
1321  {
1322  label facei = checkFaces[i];
1323 
1324  const face& f = fcs[facei];
1325 
1326  vector faceNormal = faceAreas[facei];
1327  faceNormal /= mag(faceNormal) + vSmall;
1328 
1329  // Get edge from f[0] to f[size-1];
1330  vector ePrev(p[f.first()] - p[f.last()]);
1331  scalar magEPrev = mag(ePrev);
1332  ePrev /= magEPrev + vSmall;
1333 
1334  forAll(f, fp0)
1335  {
1336  // Get vertex after fp
1337  label fp1 = f.fcIndex(fp0);
1338 
1339  // Normalised vector between two consecutive points
1340  vector e10(p[f[fp1]] - p[f[fp0]]);
1341  scalar magE10 = mag(e10);
1342  e10 /= magE10 + vSmall;
1343 
1344  if (magEPrev > small && magE10 > small)
1345  {
1346  vector edgeNormal = ePrev ^ e10;
1347  scalar magEdgeNormal = mag(edgeNormal);
1348 
1349  if (magEdgeNormal < maxSin)
1350  {
1351  // Edges (almost) aligned -> face is ok.
1352  }
1353  else
1354  {
1355  // Check normal
1356  edgeNormal /= magEdgeNormal;
1357 
1358  if ((edgeNormal & faceNormal) < small)
1359  {
1360  if (facei != errorFacei)
1361  {
1362  // Count only one error per face.
1363  errorFacei = facei;
1364  nConcave++;
1365  }
1366 
1367  if (setPtr)
1368  {
1369  setPtr->insert(facei);
1370  }
1371 
1372  maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
1373  }
1374  }
1375  }
1376 
1377  ePrev = e10;
1378  magEPrev = magE10;
1379  }
1380  }
1381 
1382  reduce(nConcave, sumOp<label>());
1383  reduce(maxEdgeSin, maxOp<scalar>());
1384 
1385  if (report)
1386  {
1387  if (maxEdgeSin > small)
1388  {
1389  scalar maxConcaveDegr =
1390  radToDeg(Foam::asin(Foam::min(1.0, maxEdgeSin)));
1391 
1392  Info<< "There are " << nConcave
1393  << " faces with concave angles between consecutive"
1394  << " edges. Max concave angle = " << maxConcaveDegr
1395  << " degrees.\n" << endl;
1396  }
1397  else
1398  {
1399  Info<< "All angles in faces are convex or less than " << maxDeg
1400  << " degrees concave.\n" << endl;
1401  }
1402  }
1403 
1404  if (nConcave > 0)
1405  {
1406  if (report)
1407  {
1409  << nConcave << " face points with severe concave angle (> "
1410  << maxDeg << " deg) found.\n"
1411  << endl;
1412  }
1413 
1414  return true;
1415  }
1416  else
1417  {
1418  return false;
1419  }
1420 }
1421 
1422 
1425  const bool report,
1426  const scalar minTwist,
1427  const polyMesh& mesh,
1428  const vectorField& cellCentres,
1429  const vectorField& faceAreas,
1430  const vectorField& faceCentres,
1431  const pointField& p,
1432  const labelList& checkFaces,
1433  labelHashSet* setPtr
1434 )
1435 {
1436  if (minTwist < -1-small || minTwist > 1+small)
1437  {
1439  << "minTwist should be [-1..1] but is now " << minTwist
1440  << abort(FatalError);
1441  }
1442 
1443 
1444  const faceList& fcs = mesh.faces();
1445 
1446 
1447  label nWarped = 0;
1448 
1449 // forAll(checkFaces, i)
1450 // {
1451 // label facei = checkFaces[i];
1452 //
1453 // const face& f = fcs[facei];
1454 //
1455 // scalar magArea = mag(faceAreas[facei]);
1456 //
1457 // if (f.size() > 3 && magArea > vSmall)
1458 // {
1459 // const vector nf = faceAreas[facei] / magArea;
1460 //
1461 // const point& fc = faceCentres[facei];
1462 //
1463 // forAll(f, fpI)
1464 // {
1465 // vector triArea
1466 // (
1467 // triPointRef
1468 // (
1469 // p[f[fpI]],
1470 // p[f.nextLabel(fpI)],
1471 // fc
1472 // ).area()
1473 // );
1474 //
1475 // scalar magTri = mag(triArea);
1476 //
1477 // if (magTri > vSmall && ((nf & triArea/magTri) < minTwist))
1478 // {
1479 // nWarped++;
1480 //
1481 // if (setPtr)
1482 // {
1483 // setPtr->insert(facei);
1484 // }
1485 //
1486 // break;
1487 // }
1488 // }
1489 // }
1490 // }
1491 
1492 
1493  const labelList& own = mesh.faceOwner();
1494  const labelList& nei = mesh.faceNeighbour();
1495  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1496 
1497  // Calculate coupled cell centre
1498  pointField neiCc(mesh.nFaces()-mesh.nInternalFaces());
1499 
1500  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
1501  {
1502  neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
1503  }
1505 
1506  forAll(checkFaces, i)
1507  {
1508  label facei = checkFaces[i];
1509 
1510  const face& f = fcs[facei];
1511 
1512  if (f.size() > 3)
1513  {
1514  vector nf(Zero);
1515 
1516  if (mesh.isInternalFace(facei))
1517  {
1518  nf = cellCentres[nei[facei]] - cellCentres[own[facei]];
1519  nf /= mag(nf) + vSmall;
1520  }
1521  else if (patches[patches.whichPatch(facei)].coupled())
1522  {
1523  nf =
1524  neiCc[facei-mesh.nInternalFaces()]
1525  - cellCentres[own[facei]];
1526  nf /= mag(nf) + vSmall;
1527  }
1528  else
1529  {
1530  nf = faceCentres[facei] - cellCentres[own[facei]];
1531  nf /= mag(nf) + vSmall;
1532  }
1533 
1534  if (nf != vector::zero)
1535  {
1536  const point& fc = faceCentres[facei];
1537 
1538  forAll(f, fpI)
1539  {
1540  vector triArea
1541  (
1542  triPointRef
1543  (
1544  p[f[fpI]],
1545  p[f.nextLabel(fpI)],
1546  fc
1547  ).area()
1548  );
1549 
1550  scalar magTri = mag(triArea);
1551 
1552  if (magTri > vSmall && ((nf & triArea/magTri) < minTwist))
1553  {
1554  nWarped++;
1555 
1556  if (setPtr)
1557  {
1558  setPtr->insert(facei);
1559  }
1560 
1561  break;
1562  }
1563  }
1564  }
1565  }
1566  }
1567 
1568  reduce(nWarped, sumOp<label>());
1569 
1570  if (report)
1571  {
1572  if (nWarped> 0)
1573  {
1574  Info<< "There are " << nWarped
1575  << " faces with cosine of the angle"
1576  << " between triangle normal and face normal less than "
1577  << minTwist << nl << endl;
1578  }
1579  else
1580  {
1581  Info<< "All faces are flat in that the cosine of the angle"
1582  << " between triangle normal and face normal less than "
1583  << minTwist << nl << endl;
1584  }
1585  }
1586 
1587  if (nWarped > 0)
1588  {
1589  if (report)
1590  {
1592  << nWarped << " faces with severe warpage "
1593  << "(cosine of the angle between triangle normal and "
1594  << "face normal < " << minTwist << ") found.\n"
1595  << endl;
1596  }
1597 
1598  return true;
1599  }
1600  else
1601  {
1602  return false;
1603  }
1604 }
1605 
1606 
1609  const bool report,
1610  const scalar minTwist,
1611  const polyMesh& mesh,
1612  const vectorField& faceAreas,
1613  const vectorField& faceCentres,
1614  const pointField& p,
1615  const labelList& checkFaces,
1616  labelHashSet* setPtr
1617 )
1618 {
1619  if (minTwist < -1-small || minTwist > 1+small)
1620  {
1622  << "minTwist should be [-1..1] but is now " << minTwist
1623  << abort(FatalError);
1624  }
1625 
1626  const faceList& fcs = mesh.faces();
1627 
1628  label nWarped = 0;
1629 
1630  forAll(checkFaces, i)
1631  {
1632  label facei = checkFaces[i];
1633 
1634  const face& f = fcs[facei];
1635 
1636  if (f.size() > 3)
1637  {
1638  const point& fc = faceCentres[facei];
1639 
1640  // Find starting triangle (at startFp) with non-zero area
1641  label startFp = -1;
1642  vector prevN;
1643 
1644  forAll(f, fp)
1645  {
1646  prevN = triPointRef
1647  (
1648  p[f[fp]],
1649  p[f.nextLabel(fp)],
1650  fc
1651  ).area();
1652 
1653  scalar magTri = mag(prevN);
1654 
1655  if (magTri > vSmall)
1656  {
1657  startFp = fp;
1658  prevN /= magTri;
1659  break;
1660  }
1661  }
1662 
1663  if (startFp != -1)
1664  {
1665  label fp = startFp;
1666 
1667  do
1668  {
1669  fp = f.fcIndex(fp);
1670 
1671  vector triN
1672  (
1673  triPointRef
1674  (
1675  p[f[fp]],
1676  p[f.nextLabel(fp)],
1677  fc
1678  ).area()
1679  );
1680  scalar magTri = mag(triN);
1681 
1682  if (magTri > vSmall)
1683  {
1684  triN /= magTri;
1685 
1686  if ((prevN & triN) < minTwist)
1687  {
1688  nWarped++;
1689 
1690  if (setPtr)
1691  {
1692  setPtr->insert(facei);
1693  }
1694 
1695  break;
1696  }
1697 
1698  prevN = triN;
1699  }
1700  else if (minTwist > 0)
1701  {
1702  nWarped++;
1703 
1704  if (setPtr)
1705  {
1706  setPtr->insert(facei);
1707  }
1708 
1709  break;
1710  }
1711  }
1712  while (fp != startFp);
1713  }
1714  }
1715  }
1716 
1717 
1718  reduce(nWarped, sumOp<label>());
1719 
1720  if (report)
1721  {
1722  if (nWarped> 0)
1723  {
1724  Info<< "There are " << nWarped
1725  << " faces with cosine of the angle"
1726  << " between consecutive triangle normals less than "
1727  << minTwist << nl << endl;
1728  }
1729  else
1730  {
1731  Info<< "All faces are flat in that the cosine of the angle"
1732  << " between consecutive triangle normals is less than "
1733  << minTwist << nl << endl;
1734  }
1735  }
1736 
1737  if (nWarped > 0)
1738  {
1739  if (report)
1740  {
1742  << nWarped << " faces with severe warpage "
1743  << "(cosine of the angle between consecutive triangle normals"
1744  << " < " << minTwist << ") found.\n"
1745  << endl;
1746  }
1747 
1748  return true;
1749  }
1750  else
1751  {
1752  return false;
1753  }
1754 }
1755 
1756 
1759  const bool report,
1760  const scalar minFlatness,
1761  const polyMesh& mesh,
1762  const vectorField& faceAreas,
1763  const vectorField& faceCentres,
1764  const pointField& p,
1765  const labelList& checkFaces,
1766  labelHashSet* setPtr
1767 )
1768 {
1769  if (minFlatness < -small || minFlatness > 1+small)
1770  {
1772  << "minFlatness should be [0..1] but is now " << minFlatness
1773  << abort(FatalError);
1774  }
1775 
1776  const faceList& fcs = mesh.faces();
1777 
1778  label nWarped = 0;
1779 
1780  forAll(checkFaces, i)
1781  {
1782  label facei = checkFaces[i];
1783 
1784  const face& f = fcs[facei];
1785 
1786  if (f.size() > 3)
1787  {
1788  const point& fc = faceCentres[facei];
1789 
1790  // Sum triangle areas
1791  scalar sumArea = 0.0;
1792 
1793  forAll(f, fp)
1794  {
1795  sumArea += triPointRef
1796  (
1797  p[f[fp]],
1798  p[f.nextLabel(fp)],
1799  fc
1800  ).mag();
1801  }
1802 
1803  if (sumArea/mag(faceAreas[facei]) < minFlatness)
1804  {
1805  nWarped++;
1806 
1807  if (setPtr)
1808  {
1809  setPtr->insert(facei);
1810  }
1811  }
1812  }
1813  }
1814 
1815  reduce(nWarped, sumOp<label>());
1816 
1817  if (report)
1818  {
1819  if (nWarped> 0)
1820  {
1821  Info<< "There are " << nWarped
1822  << " faces with area of individual triangles"
1823  << " compared to overall area less than "
1824  << minFlatness << nl << endl;
1825  }
1826  else
1827  {
1828  Info<< "All faces are flat in that the area of individual triangles"
1829  << " compared to overall area is less than "
1830  << minFlatness << nl << endl;
1831  }
1832  }
1833 
1834  if (nWarped > 0)
1835  {
1836  if (report)
1837  {
1839  << nWarped << " non-flat faces "
1840  << "(area of individual triangles"
1841  << " compared to overall area"
1842  << " < " << minFlatness << ") found.\n"
1843  << endl;
1844  }
1845 
1846  return true;
1847  }
1848  else
1849  {
1850  return false;
1851  }
1852 }
1853 
1854 
1857  const bool report,
1858  const scalar minArea,
1859  const polyMesh& mesh,
1860  const vectorField& faceAreas,
1861  const labelList& checkFaces,
1862  labelHashSet* setPtr
1863 )
1864 {
1865  label nZeroArea = 0;
1866 
1867  forAll(checkFaces, i)
1868  {
1869  label facei = checkFaces[i];
1870 
1871  if (mag(faceAreas[facei]) < minArea)
1872  {
1873  if (setPtr)
1874  {
1875  setPtr->insert(facei);
1876  }
1877  nZeroArea++;
1878  }
1879  }
1880 
1881 
1882  reduce(nZeroArea, sumOp<label>());
1883 
1884  if (report)
1885  {
1886  if (nZeroArea > 0)
1887  {
1888  Info<< "There are " << nZeroArea
1889  << " faces with area < " << minArea << '.' << nl << endl;
1890  }
1891  else
1892  {
1893  Info<< "All faces have area > " << minArea << '.' << nl << endl;
1894  }
1895  }
1896 
1897  if (nZeroArea > 0)
1898  {
1899  if (report)
1900  {
1902  << nZeroArea << " faces with area < " << minArea
1903  << " found.\n"
1904  << endl;
1905  }
1906 
1907  return true;
1908  }
1909  else
1910  {
1911  return false;
1912  }
1913 }
1914 
1915 
1918  const bool report,
1919  const scalar warnDet,
1920  const polyMesh& mesh,
1921  const vectorField& faceAreas,
1922  const labelList& checkFaces,
1923  labelHashSet* setPtr
1924 )
1925 {
1926  const cellList& cells = mesh.cells();
1927 
1928  scalar minDet = great;
1929  scalar sumDet = 0.0;
1930  label nSumDet = 0;
1931  label nWarnDet = 0;
1932 
1933  const labelList affectedCells(getAffectedCells(mesh, checkFaces));
1934 
1935  forAll(affectedCells, i)
1936  {
1937  const cell& cFaces = cells[affectedCells[i]];
1938 
1939  tensor areaSum(Zero);
1940  scalar magAreaSum = 0;
1941 
1942  forAll(cFaces, cFacei)
1943  {
1944  label facei = cFaces[cFacei];
1945 
1946  scalar magArea = mag(faceAreas[facei]);
1947 
1948  magAreaSum += magArea;
1949  areaSum += faceAreas[facei]*(faceAreas[facei]/(magArea+vSmall));
1950  }
1951 
1952  scalar scaledDet = det(areaSum/(magAreaSum+vSmall))/0.037037037037037;
1953 
1954  minDet = min(minDet, scaledDet);
1955  sumDet += scaledDet;
1956  nSumDet++;
1957 
1958  if (scaledDet < warnDet)
1959  {
1960  if (setPtr)
1961  {
1962  // Insert all faces of the cell.
1963  forAll(cFaces, cFacei)
1964  {
1965  label facei = cFaces[cFacei];
1966  setPtr->insert(facei);
1967  }
1968  }
1969  nWarnDet++;
1970  }
1971  }
1972 
1973  reduce(minDet, minOp<scalar>());
1974  reduce(sumDet, sumOp<scalar>());
1975  reduce(nSumDet, sumOp<label>());
1976  reduce(nWarnDet, sumOp<label>());
1977 
1978  if (report)
1979  {
1980  if (nSumDet > 0)
1981  {
1982  Info<< "Cell determinant (1 = uniform cube) : average = "
1983  << sumDet / nSumDet << " min = " << minDet << endl;
1984  }
1985 
1986  if (nWarnDet > 0)
1987  {
1988  Info<< "There are " << nWarnDet
1989  << " cells with determinant < " << warnDet << '.' << nl
1990  << endl;
1991  }
1992  else
1993  {
1994  Info<< "All faces have determinant > " << warnDet << '.' << nl
1995  << endl;
1996  }
1997  }
1998 
1999  if (nWarnDet > 0)
2000  {
2001  if (report)
2002  {
2004  << nWarnDet << " cells with determinant < " << warnDet
2005  << " found.\n"
2006  << endl;
2007  }
2008 
2009  return true;
2010  }
2011  else
2012  {
2013  return false;
2014  }
2015 }
2016 
2017 
2018 // ************************************************************************* //
const fvPatchList & patches
bool checkTriangleTwist(const bool report, const scalar minTwist, const polyMesh &, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Check consecutive face-triangle (from face-centre decomposition) normals.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
dimensionedScalar acos(const dimensionedScalar &ds)
tetrahedron< point, const point & > tetPointRef
Definition: tetPointRef.H:44
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A triangle primitive used to calculate face areas and swept volumes.
Definition: triangle.H:57
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
bool checkFaceArea(const bool report, const scalar minArea, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
Check for small faces.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
scalar checkNonOrtho(const polyMesh &mesh, const bool report, const scalar severeNonorthogonalityThreshold, const label facei, const vector &s, const vector &d, label &severeNonOrth, label &errorNonOrth, labelHashSet *setPtr)
Definition: polyMeshCheck.C:44
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
label nInternalFaces() const
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1243
static void swapBoundaryFacePositions(const polyMesh &mesh, UList< point > &l)
Swap coupled positions.
Definition: syncTools.H:452
bool checkFaceTet(const polyMesh &mesh, const bool report, const scalar minTetQuality, const pointField &p, const label facei, const point &fc, const point &cc, labelHashSet *setPtr)
label nFaces() const
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:100
const cellList & cells() const
dimensionedScalar det(const dimensionedSphericalTensor &dt)
static void swapBoundaryCellPositions(const polyMesh &mesh, const UList< point > &cellData, List< point > &neighbourCellData)
Swap to obtain neighbour cell positions for all boundary faces.
Definition: syncTools.C:31
bool checkFaceDotProduct(const bool report, const scalar orthWarn, const polyMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
Check face non-orthogonality.
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
T & first()
Return the first element of the list.
Definition: UListI.H:114
bool checkFaceSkewness(const bool report, const scalar internalSkew, const scalar boundarySkew, const polyMesh &mesh, const pointField &points, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
Check face skewness.
dimensionedScalar asin(const dimensionedScalar &ds)
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
static label findBasePoint(const polyMesh &mesh, label fI, scalar tol, bool report=false)
Find the base point to use for a minimum triangle.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
bool checkFaceTets(const bool report, const scalar minPyrVol, const polyMesh &, const vectorField &cellCentres, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *)
Check face tetrahedra volumes.
static tmp< scalarField > faceSkewness(const primitiveMesh &mesh, const pointField &points, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs)
Generate skewness field.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
static scalar boundaryFaceSkewness(const primitiveMesh &mesh, const pointField &p, const vectorField &fCtrs, const vectorField &fAreas, const label facei, const point &ownCc)
Skewness of single boundary face.
dimensionedScalar cos(const dimensionedScalar &ds)
const cellShapeList & cells
PolyMesh checking routines. Checks various criteria for a mesh and supplied geometry, with the mesh only used for topology. Coupled patch aware (i.e., coupled faces are treated as internal).
static label findSharedBasePoint(const polyMesh &mesh, label fI, const point &nCc, scalar tol, bool report=false)
Find the first suitable base point to use for a minimum.
labelList getAffectedCells(const polyMesh &mesh, const labelList &changedFaces)
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1237
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
static const zero Zero
Definition: zero.H:97
bool checkFaceTwist(const bool report, const scalar minTwist, const polyMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1224
pyramid< point, const point &, const face & > pyramidPointFaceRef
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::polyBoundaryMesh.
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
bool checkFaceFlatness(const bool report, const scalar minFlatness, const polyMesh &, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Check for face areas v.s. sum of face-triangle (from face-centre.
dimensionedScalar sin(const dimensionedScalar &ds)
bool checkFaceAngles(const bool report, const scalar maxDeg, const polyMesh &mesh, const vectorField &faceAreas, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Check convexity of angles in a face. See primitiveMesh for explanation.
static const char nl
Definition: Ostream.H:260
labelList f(nPoints)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
bool checkCellDeterminant(const bool report, const scalar minDet, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
Check the area of internal faces v.s. boundary faces.
label patchi
bool checkFacePyramids(const bool report, const scalar minPyrVol, const polyMesh &, const vectorField &cellCentres, const pointField &p, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *)
Check face pyramid volumes.
#define WarningInFunction
Report a warning using Foam::Warning.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
bool checkFaceWeights(const bool report, const scalar warnWeight, const polyMesh &mesh, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
Check interpolation weights (0.5 for regular mesh)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:436
T & last()
Return the last element of the list.
Definition: UListI.H:128
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
bool checkVolRatio(const bool report, const scalar warnRatio, const polyMesh &mesh, const scalarField &cellVolumes, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
Cell volume ratio of neighbouring cells (1 for regular mesh)
Namespace for OpenFOAM.