motionSmootherAlgoCheck.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2014 OpenFOAM Foundation
9  Copyright (C) 2020,2022 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "motionSmootherAlgo.H"
30 #include "polyMeshGeometry.H"
31 #include "IOmanip.H"
32 #include "pointSet.H"
33 
34 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
35 
36 bool Foam::motionSmootherAlgo::checkMesh
37 (
38  const bool report,
39  const polyMesh& mesh,
40  const dictionary& dict,
41  const labelList& checkFaces,
42  labelHashSet& wrongFaces,
43  const bool dryRun
44 )
45 {
46  List<labelPair> emptyBaffles;
47  return checkMesh
48  (
49  report,
50  mesh,
51  dict,
52  checkFaces,
53  emptyBaffles,
54  wrongFaces,
55  dryRun
56  );
57 }
58 
59 bool Foam::motionSmootherAlgo::checkMesh
60 (
61  const bool report,
62  const polyMesh& mesh,
63  const dictionary& dict,
64  const labelList& checkFaces,
65  const List<labelPair>& baffles,
66  labelHashSet& wrongFaces,
67  const bool dryRun
68 )
69 {
70  const scalar maxNonOrtho
71  (
72  get<scalar>(dict, "maxNonOrtho", dryRun, keyType::REGEX_RECURSIVE)
73  );
74  const scalar minVol
75  (
76  get<scalar>(dict, "minVol", dryRun, keyType::REGEX_RECURSIVE)
77  );
78  const scalar minTetQuality
79  (
80  get<scalar>(dict, "minTetQuality", dryRun, keyType::REGEX_RECURSIVE)
81  );
82  const scalar maxConcave
83  (
84  get<scalar>(dict, "maxConcave", dryRun, keyType::REGEX_RECURSIVE)
85  );
86  const scalar minArea
87  (
88  get<scalar>(dict, "minArea", dryRun, keyType::REGEX_RECURSIVE)
89  );
90  const scalar maxIntSkew
91  (
92  get<scalar>
93  (
94  dict, "maxInternalSkewness", dryRun, keyType::REGEX_RECURSIVE
95  )
96  );
97  const scalar maxBounSkew
98  (
99  get<scalar>
100  (
101  dict, "maxBoundarySkewness", dryRun, keyType::REGEX_RECURSIVE
102  )
103  );
104  const scalar minWeight
105  (
106  get<scalar>(dict, "minFaceWeight", dryRun, keyType::REGEX_RECURSIVE)
107  );
108  const scalar minVolRatio
109  (
110  get<scalar>(dict, "minVolRatio", dryRun, keyType::REGEX_RECURSIVE)
111  );
112  const scalar minTwist
113  (
114  get<scalar>(dict, "minTwist", dryRun, keyType::REGEX_RECURSIVE)
115  );
116  const scalar minTriangleTwist
117  (
118  get<scalar>(dict, "minTriangleTwist", dryRun, keyType::REGEX_RECURSIVE)
119  );
120  const scalar minFaceFlatness
121  (
122  dict.getOrDefault<scalar>
123  (
124  "minFaceFlatness", -1, keyType::REGEX_RECURSIVE
125  )
126  );
127  const scalar minDet
128  (
129  get<scalar>(dict, "minDeterminant", dryRun, keyType::REGEX_RECURSIVE)
130  );
131  const scalar minEdgeLength
132  (
133  dict.getOrDefault<scalar>
134  (
135  "minEdgeLength", -1, keyType::REGEX_RECURSIVE
136  )
137  );
138 
139 
140  if (dryRun)
141  {
142  string errorMsg(FatalError.message());
143  string IOerrorMsg(FatalIOError.message());
144 
145  if (errorMsg.size() || IOerrorMsg.size())
146  {
147  //errorMsg = "[dryRun] " + errorMsg;
148  //errorMsg.replaceAll("\n", "\n[dryRun] ");
149  //IOerrorMsg = "[dryRun] " + IOerrorMsg;
150  //IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
151 
152  IOWarningInFunction(dict)
153  << nl
154  << "Missing/incorrect required dictionary entries:" << nl
155  << nl
156  << IOerrorMsg.c_str() << nl
157  << errorMsg.c_str() << nl
158  //<< nl << "Exiting dry-run" << nl
159  << endl;
160 
161  FatalError.clear();
162  FatalIOError.clear();
163  }
164  return false;
165  }
166 
167 
168  label nWrongFaces = 0;
169 
170  Info<< "Checking faces in error :" << endl;
171  //Pout.setf(ios_base::left);
172 
173  if (maxNonOrtho < 180.0-SMALL)
174  {
175  polyMeshGeometry::checkFaceDotProduct
176  (
177  report,
178  maxNonOrtho,
179  mesh,
180  mesh.cellCentres(),
181  mesh.faceAreas(),
182  checkFaces,
183  baffles,
184  &wrongFaces
185  );
186 
187  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
188 
189  Info<< " non-orthogonality > "
190  << setw(3) << maxNonOrtho
191  << " degrees : "
192  << nNewWrongFaces-nWrongFaces << endl;
193 
194  nWrongFaces = nNewWrongFaces;
195  }
196 
197  if (minVol > -GREAT)
198  {
199  polyMeshGeometry::checkFacePyramids
200  (
201  report,
202  minVol,
203  mesh,
204  mesh.cellCentres(),
205  mesh.points(),
206  checkFaces,
207  baffles,
208  &wrongFaces
209  );
210 
211  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
212 
213  Info<< " faces with face pyramid volume < "
214  << setw(5) << minVol << " : "
215  << nNewWrongFaces-nWrongFaces << endl;
216 
217  nWrongFaces = nNewWrongFaces;
218  }
219 
220  if (minTetQuality > -GREAT)
221  {
222  polyMeshGeometry::checkFaceTets
223  (
224  report,
225  minTetQuality,
226  mesh,
227  mesh.cellCentres(),
228  mesh.faceCentres(),
229  mesh.points(),
230  checkFaces,
231  baffles,
232  &wrongFaces
233  );
234 
235  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
236 
237  Info<< " faces with face-decomposition tet quality < "
238  << setw(5) << minTetQuality << " : "
239  << nNewWrongFaces-nWrongFaces << endl;
240 
241  nWrongFaces = nNewWrongFaces;
242  }
243 
244  if (maxConcave < 180.0-SMALL)
245  {
246  polyMeshGeometry::checkFaceAngles
247  (
248  report,
249  maxConcave,
250  mesh,
251  mesh.faceAreas(),
252  mesh.points(),
253  checkFaces,
254  &wrongFaces
255  );
256 
257  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
258 
259  Info<< " faces with concavity > "
260  << setw(3) << maxConcave
261  << " degrees : "
262  << nNewWrongFaces-nWrongFaces << endl;
263 
264  nWrongFaces = nNewWrongFaces;
265  }
266 
267  if (minArea > -SMALL)
268  {
269  polyMeshGeometry::checkFaceArea
270  (
271  report,
272  minArea,
273  mesh,
274  mesh.faceAreas(),
275  checkFaces,
276  &wrongFaces
277  );
278 
279  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
280 
281  Info<< " faces with area < "
282  << setw(5) << minArea
283  << " m^2 : "
284  << nNewWrongFaces-nWrongFaces << endl;
285 
286  nWrongFaces = nNewWrongFaces;
287  }
288 
289  if (maxIntSkew > 0 || maxBounSkew > 0)
290  {
291  polyMeshGeometry::checkFaceSkewness
292  (
293  report,
294  maxIntSkew,
295  maxBounSkew,
296  mesh,
297  mesh.points(),
298  mesh.cellCentres(),
299  mesh.faceCentres(),
300  mesh.faceAreas(),
301  checkFaces,
302  baffles,
303  &wrongFaces
304  );
305 
306  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
307 
308  Info<< " faces with skewness > "
309  << setw(3) << maxIntSkew
310  << " (internal) or " << setw(3) << maxBounSkew
311  << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
312 
313  nWrongFaces = nNewWrongFaces;
314  }
315 
316  if (minWeight >= 0 && minWeight < 1)
317  {
318  polyMeshGeometry::checkFaceWeights
319  (
320  report,
321  minWeight,
322  mesh,
323  mesh.cellCentres(),
324  mesh.faceCentres(),
325  mesh.faceAreas(),
326  checkFaces,
327  baffles,
328  &wrongFaces
329  );
330 
331  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
332 
333  Info<< " faces with interpolation weights (0..1) < "
334  << setw(5) << minWeight
335  << " : "
336  << nNewWrongFaces-nWrongFaces << endl;
337 
338  nWrongFaces = nNewWrongFaces;
339  }
340 
341  if (minVolRatio >= 0)
342  {
343  polyMeshGeometry::checkVolRatio
344  (
345  report,
346  minVolRatio,
347  mesh,
348  mesh.cellVolumes(),
349  checkFaces,
350  baffles,
351  &wrongFaces
352  );
353 
354  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
355 
356  Info<< " faces with volume ratio of neighbour cells < "
357  << setw(5) << minVolRatio
358  << " : "
359  << nNewWrongFaces-nWrongFaces << endl;
360 
361  nWrongFaces = nNewWrongFaces;
362  }
363 
364  if (minTwist > -1)
365  {
366  //Pout<< "Checking face twist: dot product of face normal "
367  // << "with face triangle normals" << endl;
368  polyMeshGeometry::checkFaceTwist
369  (
370  report,
371  minTwist,
372  mesh,
373  mesh.cellCentres(),
374  mesh.faceAreas(),
375  mesh.faceCentres(),
376  mesh.points(),
377  checkFaces,
378  &wrongFaces
379  );
380 
381  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
382 
383  Info<< " faces with face twist < "
384  << setw(5) << minTwist
385  << " : "
386  << nNewWrongFaces-nWrongFaces << endl;
387 
388  nWrongFaces = nNewWrongFaces;
389  }
390 
391  if (minTriangleTwist > -1)
392  {
393  //Pout<< "Checking triangle twist: dot product of consecutive triangle"
394  // << " normals resulting from face-centre decomposition" << endl;
395  polyMeshGeometry::checkTriangleTwist
396  (
397  report,
398  minTriangleTwist,
399  mesh,
400  mesh.faceAreas(),
401  mesh.faceCentres(),
402  mesh.points(),
403  checkFaces,
404  &wrongFaces
405  );
406 
407  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
408 
409  Info<< " faces with triangle twist < "
410  << setw(5) << minTriangleTwist
411  << " : "
412  << nNewWrongFaces-nWrongFaces << endl;
413 
414  nWrongFaces = nNewWrongFaces;
415  }
416 
417  if (minFaceFlatness > -SMALL)
418  {
419  polyMeshGeometry::checkFaceFlatness
420  (
421  report,
422  minFaceFlatness,
423  mesh,
424  mesh.faceAreas(),
425  mesh.faceCentres(),
426  mesh.points(),
427  checkFaces,
428  &wrongFaces
429  );
430 
431  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
432 
433  Info<< " faces with flatness < "
434  << setw(5) << minFaceFlatness
435  << " : "
436  << nNewWrongFaces-nWrongFaces << endl;
437 
438  nWrongFaces = nNewWrongFaces;
439  }
440 
441  if (minDet > -1)
442  {
443  polyMeshGeometry::checkCellDeterminant
444  (
445  report,
446  minDet,
447  mesh,
448  mesh.faceAreas(),
449  checkFaces,
450  polyMeshGeometry::affectedCells(mesh, checkFaces),
451  &wrongFaces
452  );
453 
454  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
455 
456  Info<< " faces on cells with determinant < "
457  << setw(5) << minDet << " : "
458  << nNewWrongFaces-nWrongFaces << endl;
459 
460  nWrongFaces = nNewWrongFaces;
461  }
462 
463  if (minEdgeLength >= 0)
464  {
465  pointSet edgePoints(mesh, "smallEdgePoints", mesh.nPoints()/1000);
466  polyMeshGeometry::checkEdgeLength
467  (
468  report,
469  minEdgeLength,
470  mesh,
471  checkFaces,
472  &edgePoints
473  );
474 
475  const auto& pointFaces = mesh.pointFaces();
476 
477  label nNewWrongFaces = 0;
478  for (const label pointi : edgePoints)
479  {
480  const auto& pFaces = pointFaces[pointi];
481  for (const label facei : pFaces)
482  {
483  if (wrongFaces.insert(facei))
484  {
485  nNewWrongFaces++;
486  }
487  }
488  }
489 
490  Info<< " faces with edge length < "
491  << setw(5) << minEdgeLength << " : "
492  << returnReduce(nNewWrongFaces, sumOp<label>()) << endl;
493 
494  nWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
495  }
496 
497  //Pout.setf(ios_base::right);
498 
499  return nWrongFaces > 0;
500 }
501 
502 
503 bool Foam::motionSmootherAlgo::checkMesh
504 (
505  const bool report,
506  const polyMesh& mesh,
507  const dictionary& dict,
508  labelHashSet& wrongFaces,
509  const bool dryRun
510 )
511 {
512  return checkMesh
513  (
514  report,
515  mesh,
516  dict,
517  identity(mesh.nFaces()),
518  wrongFaces,
519  dryRun
520  );
521 }
522 
523 bool Foam::motionSmootherAlgo::checkMesh
524 (
525  const bool report,
526  const dictionary& dict,
527  const polyMeshGeometry& meshGeom,
528  const pointField& points,
529  const labelList& checkFaces,
530  labelHashSet& wrongFaces,
531  const bool dryRun
532 )
533 {
534  List<labelPair> emptyBaffles;
535 
536  return checkMesh
537  (
538  report,
539  dict,
540  meshGeom,
541  points,
542  checkFaces,
543  emptyBaffles,
544  wrongFaces,
545  dryRun
546  );
547 }
548 
549 
550 bool Foam::motionSmootherAlgo::checkMesh
551 (
552  const bool report,
553  const dictionary& dict,
554  const polyMeshGeometry& meshGeom,
555  const pointField& points,
556  const labelList& checkFaces,
557  const List<labelPair>& baffles,
558  labelHashSet& wrongFaces,
559  const bool dryRun
560 )
561 {
562  const scalar maxNonOrtho
563  (
564  get<scalar>(dict, "maxNonOrtho", dryRun, keyType::REGEX_RECURSIVE)
565  );
566  const scalar minVol
567  (
568  get<scalar>(dict, "minVol", dryRun, keyType::REGEX_RECURSIVE)
569  );
570  const scalar minTetQuality
571  (
572  get<scalar>(dict, "minTetQuality", dryRun, keyType::REGEX_RECURSIVE)
573  );
574  const scalar maxConcave
575  (
576  get<scalar>(dict, "maxConcave", dryRun, keyType::REGEX_RECURSIVE)
577  );
578  const scalar minArea
579  (
580  get<scalar>(dict, "minArea", dryRun, keyType::REGEX_RECURSIVE)
581  );
582  const scalar maxIntSkew
583  (
584  get<scalar>
585  (
586  dict,
587  "maxInternalSkewness",
588  dryRun,
589  keyType::REGEX_RECURSIVE
590  )
591  );
592  const scalar maxBounSkew
593  (
594  get<scalar>
595  (
596  dict,
597  "maxBoundarySkewness",
598  dryRun,
599  keyType::REGEX_RECURSIVE
600  )
601  );
602  const scalar minWeight
603  (
604  get<scalar>(dict, "minFaceWeight", dryRun, keyType::REGEX_RECURSIVE)
605  );
606  const scalar minVolRatio
607  (
608  get<scalar>(dict, "minVolRatio", dryRun, keyType::REGEX_RECURSIVE)
609  );
610  const scalar minTwist
611  (
612  get<scalar>(dict, "minTwist", dryRun, keyType::REGEX_RECURSIVE)
613  );
614  const scalar minTriangleTwist
615  (
616  get<scalar>(dict, "minTriangleTwist", dryRun, keyType::REGEX_RECURSIVE)
617  );
618  scalar minFaceFlatness = -1.0;
619  dict.readIfPresent
620  (
621  "minFaceFlatness",
622  minFaceFlatness,
623  keyType::REGEX_RECURSIVE
624  );
625  const scalar minDet
626  (
627  get<scalar>(dict, "minDeterminant", dryRun, keyType::REGEX_RECURSIVE)
628  );
629  const scalar minEdgeLength
630  (
631  dict.getOrDefault<scalar>
632  (
633  "minEdgeLength", -1, keyType::REGEX_RECURSIVE
634  )
635  );
636 
637  if (dryRun)
638  {
639  string errorMsg(FatalError.message());
640  string IOerrorMsg(FatalIOError.message());
641 
642  if (errorMsg.size() || IOerrorMsg.size())
643  {
644  //errorMsg = "[dryRun] " + errorMsg;
645  //errorMsg.replaceAll("\n", "\n[dryRun] ");
646  //IOerrorMsg = "[dryRun] " + IOerrorMsg;
647  //IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
648 
649  Perr<< nl
650  << "Missing/incorrect required dictionary entries:" << nl
651  << nl
652  << IOerrorMsg.c_str() << nl
653  << errorMsg.c_str() << nl
654  //<< nl << "Exiting dry-run" << nl
655  << endl;
656 
657  FatalError.clear();
658  FatalIOError.clear();
659  }
660  return false;
661  }
662 
663 
664  label nWrongFaces = 0;
665 
666  Info<< "Checking faces in error :" << endl;
667  //Pout.setf(ios_base::left);
668 
669  if (maxNonOrtho < 180.0-SMALL)
670  {
671  meshGeom.checkFaceDotProduct
672  (
673  report,
674  maxNonOrtho,
675  checkFaces,
676  baffles,
677  &wrongFaces
678  );
679 
680  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
681 
682  Info<< " non-orthogonality > "
683  << setw(3) << maxNonOrtho
684  << " degrees : "
685  << nNewWrongFaces-nWrongFaces << endl;
686 
687  nWrongFaces = nNewWrongFaces;
688  }
689 
690  if (minVol > -GREAT)
691  {
692  meshGeom.checkFacePyramids
693  (
694  report,
695  minTetQuality,
696  points,
697  checkFaces,
698  baffles,
699  &wrongFaces
700  );
701 
702  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
703 
704  Info<< " faces with face pyramid volume < "
705  << setw(5) << minVol << " : "
706  << nNewWrongFaces-nWrongFaces << endl;
707 
708  nWrongFaces = nNewWrongFaces;
709  }
710 
711  if (minTetQuality > -GREAT)
712  {
713  meshGeom.checkFaceTets
714  (
715  report,
716  minTetQuality,
717  points,
718  checkFaces,
719  baffles,
720  &wrongFaces
721  );
722 
723  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
724 
725  Info<< " faces with face-decomposition tet quality < "
726  << setw(5) << minTetQuality << " : "
727  << nNewWrongFaces-nWrongFaces << endl;
728 
729  nWrongFaces = nNewWrongFaces;
730  }
731 
732  if (maxConcave < 180.0-SMALL)
733  {
734  meshGeom.checkFaceAngles
735  (
736  report,
737  maxConcave,
738  points,
739  checkFaces,
740  &wrongFaces
741  );
742 
743  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
744 
745  Info<< " faces with concavity > "
746  << setw(3) << maxConcave
747  << " degrees : "
748  << nNewWrongFaces-nWrongFaces << endl;
749 
750  nWrongFaces = nNewWrongFaces;
751  }
752 
753  if (minArea > -SMALL)
754  {
755  meshGeom.checkFaceArea
756  (
757  report,
758  minArea,
759  checkFaces,
760  &wrongFaces
761  );
762 
763  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
764 
765  Info<< " faces with area < "
766  << setw(5) << minArea
767  << " m^2 : "
768  << nNewWrongFaces-nWrongFaces << endl;
769 
770  nWrongFaces = nNewWrongFaces;
771  }
772 
773  if (maxIntSkew > 0 || maxBounSkew > 0)
774  {
775  polyMeshGeometry::checkFaceSkewness
776  (
777  report,
778  maxIntSkew,
779  maxBounSkew,
780  meshGeom.mesh(),
781  points,
782  meshGeom.cellCentres(),
783  meshGeom.faceCentres(),
784  meshGeom.faceAreas(),
785  checkFaces,
786  baffles,
787  &wrongFaces
788  );
789 
790  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
791 
792  Info<< " faces with skewness > "
793  << setw(3) << maxIntSkew
794  << " (internal) or " << setw(3) << maxBounSkew
795  << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
796 
797  nWrongFaces = nNewWrongFaces;
798  }
799 
800  if (minWeight >= 0 && minWeight < 1)
801  {
802  meshGeom.checkFaceWeights
803  (
804  report,
805  minWeight,
806  checkFaces,
807  baffles,
808  &wrongFaces
809  );
810 
811  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
812 
813  Info<< " faces with interpolation weights (0..1) < "
814  << setw(5) << minWeight
815  << " : "
816  << nNewWrongFaces-nWrongFaces << endl;
817 
818  nWrongFaces = nNewWrongFaces;
819  }
820 
821  if (minVolRatio >= 0)
822  {
823  meshGeom.checkVolRatio
824  (
825  report,
826  minVolRatio,
827  checkFaces,
828  baffles,
829  &wrongFaces
830  );
831 
832  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
833 
834  Info<< " faces with volume ratio of neighbour cells < "
835  << setw(5) << minVolRatio
836  << " : "
837  << nNewWrongFaces-nWrongFaces << endl;
838 
839  nWrongFaces = nNewWrongFaces;
840  }
841 
842  if (minTwist > -1)
843  {
844  //Pout<< "Checking face twist: dot product of face normal "
845  // << "with face triangle normals" << endl;
846  meshGeom.checkFaceTwist
847  (
848  report,
849  minTwist,
850  points,
851  checkFaces,
852  &wrongFaces
853  );
854 
855  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
856 
857  Info<< " faces with face twist < "
858  << setw(5) << minTwist
859  << " : "
860  << nNewWrongFaces-nWrongFaces << endl;
861 
862  nWrongFaces = nNewWrongFaces;
863  }
864 
865  if (minTriangleTwist > -1)
866  {
867  //Pout<< "Checking triangle twist: dot product of consecutive triangle"
868  // << " normals resulting from face-centre decomposition" << endl;
869  meshGeom.checkTriangleTwist
870  (
871  report,
872  minTriangleTwist,
873  points,
874  checkFaces,
875  &wrongFaces
876  );
877 
878  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
879 
880  Info<< " faces with triangle twist < "
881  << setw(5) << minTriangleTwist
882  << " : "
883  << nNewWrongFaces-nWrongFaces << endl;
884 
885  nWrongFaces = nNewWrongFaces;
886  }
887 
888  if (minFaceFlatness > -SMALL)
889  {
890  meshGeom.checkFaceFlatness
891  (
892  report,
893  minFaceFlatness,
894  points,
895  checkFaces,
896  &wrongFaces
897  );
898 
899  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
900 
901  Info<< " faces with flatness < "
902  << setw(5) << minFaceFlatness
903  << " : "
904  << nNewWrongFaces-nWrongFaces << endl;
905 
906  nWrongFaces = nNewWrongFaces;
907  }
908 
909  if (minDet > -1)
910  {
911  meshGeom.checkCellDeterminant
912  (
913  report,
914  minDet,
915  checkFaces,
916  polyMeshGeometry::affectedCells(meshGeom.mesh(), checkFaces),
917  &wrongFaces
918  );
919 
920  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
921 
922  Info<< " faces on cells with determinant < "
923  << setw(5) << minDet << " : "
924  << nNewWrongFaces-nWrongFaces << endl;
925 
926  nWrongFaces = nNewWrongFaces;
927  }
928 
929  if (minEdgeLength >= 0)
930  {
931  pointSet edgePoints
932  (
933  meshGeom.mesh(),
934  "smallEdgePoints",
935  meshGeom.mesh().nPoints()/1000
936  );
937  meshGeom.checkEdgeLength
938  (
939  report,
940  minEdgeLength,
941  checkFaces,
942  &edgePoints
943  );
944 
945  const auto& pointFaces = meshGeom.mesh().pointFaces();
946 
947  label nNewWrongFaces = 0;
948  for (const label pointi : edgePoints)
949  {
950  const auto& pFaces = pointFaces[pointi];
951  for (const label facei : pFaces)
952  {
953  if (wrongFaces.insert(facei))
954  {
955  nNewWrongFaces++;
956  }
957  }
958  }
959 
960  Info<< " faces with edge length < "
961  << setw(5) << minEdgeLength << " : "
962  << returnReduce(nNewWrongFaces, sumOp<label>()) << endl;
963 
964  nWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
965  }
966 
967  //Pout.setf(ios_base::right);
968 
969  return nWrongFaces > 0;
970 }
971 
972 
973 // ************************************************************************* //
Foam::polyMeshGeometry::checkFacePyramids
static bool checkFacePyramids(const bool report, const scalar minPyrVol, const polyMesh &, const vectorField &cellCentres, const pointField &p, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *)
See primitiveMesh.
Definition: polyMeshGeometry.C:541
Foam::motionSmootherAlgo::checkMesh
static bool checkMesh(const bool report, const polyMesh &mesh, const dictionary &dict, labelHashSet &wrongFaces, const bool dryRun=false)
Check mesh with mesh settings in dict. Collects incorrect faces.
Definition: motionSmootherAlgoCheck.C:497
Foam::polyMeshGeometry::checkFaceTets
static bool checkFaceTets(const bool report, const scalar minPyrVol, const polyMesh &, const vectorField &cellCentres, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *)
See primitiveMesh.
Definition: polyMeshGeometry.C:712
Foam::Perr
prefixOSstream Perr
OSstream wrapped stderr (std::cerr) with parallel prefix.
dict
dictionary dict
Definition: searchingEngine.H:11
Foam::motionSmootherAlgo::mesh
const polyMesh & mesh() const
Reference to mesh.
Definition: motionSmootherAlgo.C:354
Foam::error::clear
void clear() const
Clear any messages.
Definition: error.C:337
Foam::polyMeshGeometry::checkTriangleTwist
static bool checkTriangleTwist(const bool report, const scalar minTwist, const polyMesh &, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Consecutive triangle (from face-centre decomposition) normals.
Definition: polyMeshGeometry.C:1694
Foam::polyMeshGeometry::faceAreas
const vectorField & faceAreas() const
Definition: polyMeshGeometry.H:165
Foam::primitiveMesh::nPoints
label nPoints() const noexcept
Number of mesh points.
Definition: primitiveMeshI.H:30
Foam::FatalError
error FatalError
Error stream (stdout output on all processes), with additional &#39;FOAM FATAL ERROR&#39; header text and sta...
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg...
Definition: dictionary.H:129
Foam::polyMeshGeometry::checkFaceTwist
static bool checkFaceTwist(const bool report, const scalar minTwist, const polyMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Triangle (from face-centre decomposition) normal v.s.
Definition: polyMeshGeometry.C:1548
Foam::polyMeshGeometry::cellCentres
const vectorField & cellCentres() const
Definition: polyMeshGeometry.H:173
Foam::nl
constexpr char nl
The newline &#39;\n&#39; character (0x0a)
Definition: Ostream.H:50
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:531
Foam::primitiveMesh::nFaces
label nFaces() const noexcept
Number of mesh faces.
Definition: primitiveMeshI.H:83
Foam::HashSet::insert
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:232
Foam::returnReduce
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
Definition: PstreamReduceOps.H:463
Foam::polyMeshGeometry::checkVolRatio
static bool checkVolRatio(const bool report, const scalar warnRatio, const polyMesh &mesh, const scalarField &cellVolumes, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
Cell volume ratio of neighbouring cells (1 for regular mesh)
Definition: polyMeshGeometry.C:1287
Foam::polyMeshGeometry::checkFaceFlatness
static bool checkFaceFlatness(const bool report, const scalar minFlatness, const polyMesh &, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Area of faces v.s. sum of triangle areas.
Definition: polyMeshGeometry.C:1840
Foam::polyMesh::points
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1078
Foam::HashTable::size
label size() const noexcept
The number of elements in table.
Definition: HashTable.H:342
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::identity
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
Definition: labelLists.C:44
Foam::polyMeshGeometry::checkEdgeLength
static bool checkEdgeLength(const bool report, const scalar minEdgeLength, const polyMesh &mesh, const labelList &checkFaces, labelHashSet *pointSetPtr)
Small edges. Optionally collects points of small edges.
Definition: polyMeshGeometry.C:2082
Foam::error::message
string message() const
The accumulated error message.
Definition: error.C:331
Foam::polyMeshGeometry
Updateable mesh geometry and checking routines.
Definition: polyMeshGeometry.H:51
Foam::primitiveMesh::cellCentres
const vectorField & cellCentres() const
Definition: primitiveMeshCellCentresAndVols.C:78
IOmanip.H
Istream and Ostream manipulators taking arguments.
Foam::polyMeshGeometry::checkFaceArea
static bool checkFaceArea(const bool report, const scalar minArea, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
Small faces.
Definition: polyMeshGeometry.C:1933
Foam::dictionary::readIfPresent
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Find an entry if present, and assign to T val. FatalIOError if it is found and the number of tokens i...
Definition: dictionaryTemplates.C:413
pFaces
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
Foam::primitiveMesh::faceCentres
const vectorField & faceCentres() const
Definition: primitiveMeshFaceCentresAndAreas.C:71
Foam::polyMeshGeometry::faceCentres
const vectorField & faceCentres() const
Definition: polyMeshGeometry.H:169
Foam::polyMeshGeometry::checkFaceAngles
static bool checkFaceAngles(const bool report, const scalar maxDeg, const polyMesh &mesh, const vectorField &faceAreas, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
See primitiveMesh.
Definition: polyMeshGeometry.C:1423
Foam::polyMeshGeometry::checkFaceDotProduct
static bool checkFaceDotProduct(const bool report, const scalar orthWarn, const polyMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
See primitiveMesh.
Definition: polyMeshGeometry.C:363
Foam::primitiveMesh::faceAreas
const vectorField & faceAreas() const
Definition: primitiveMeshFaceCentresAndAreas.C:83
Foam::polyMeshGeometry::affectedCells
static labelList affectedCells(const polyMesh &, const labelList &changedFaces)
Helper function: get affected cells from faces.
Definition: polyMeshGeometry.C:186
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::primitiveMesh::pointFaces
const labelListList & pointFaces() const
Definition: primitiveMeshPointFaces.C:26
Foam::polyMeshGeometry::checkFaceWeights
static bool checkFaceWeights(const bool report, const scalar warnWeight, const polyMesh &mesh, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
Interpolation weights (0.5 for regular mesh)
Definition: polyMeshGeometry.C:1139
IOWarningInFunction
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
Definition: messageStream.H:440
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Foam::setw
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
Foam::polyMeshGeometry::checkFaceSkewness
static bool checkFaceSkewness(const bool report, const scalar internalSkew, const scalar boundarySkew, const polyMesh &mesh, const pointField &points, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
See primitiveMesh.
Definition: polyMeshGeometry.C:933
Foam::dictionary::getOrDefault
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T, or return the given default value. FatalIOError if it is found and the number of...
Definition: dictionaryTemplates.C:141
Foam::polyMeshGeometry::checkCellDeterminant
static bool checkCellDeterminant(const bool report, const scalar minDet, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, const labelList &affectedCells, labelHashSet *setPtr)
Area of internal faces v.s. boundary faces.
Definition: polyMeshGeometry.C:1990
Foam::polyMeshGeometry::mesh
const polyMesh & mesh() const
Definition: polyMeshGeometry.H:160
Foam::primitiveMesh::cellVolumes
const scalarField & cellVolumes() const
Definition: primitiveMeshCellCentresAndVols.C:90
Foam::FatalIOError
IOerror FatalIOError
Error stream (stdout output on all processes), with additional &#39;FOAM FATAL IO ERROR&#39; header text and ...