FreeFOAM The Cross-Platform CFD Toolkit
motionSmootherCheck.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) 1991-2010 OpenCFD Ltd.
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 "motionSmoother.H"
28 #include <OpenFOAM/IOmanip.H>
29 
30 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
31 
33 (
34  const bool report,
35  const polyMesh& mesh,
36  const dictionary& dict,
37  const labelList& checkFaces,
38  labelHashSet& wrongFaces
39 )
40 {
41  List<labelPair> emptyBaffles;
42  return checkMesh
43  (
44  report,
45  mesh,
46  dict,
47  checkFaces,
48  emptyBaffles,
49  wrongFaces
50  );
51 }
52 
54 (
55  const bool report,
56  const polyMesh& mesh,
57  const dictionary& dict,
58  const labelList& checkFaces,
59  const List<labelPair>& baffles,
60  labelHashSet& wrongFaces
61 )
62 {
63  const scalar maxNonOrtho
64  (
65  readScalar(dict.lookup("maxNonOrtho", true))
66  );
67  const scalar minVol
68  (
69  readScalar(dict.lookup("minVol", true))
70  );
71  const scalar maxConcave
72  (
73  readScalar(dict.lookup("maxConcave", true))
74  );
75  const scalar minArea
76  (
77  readScalar(dict.lookup("minArea", true))
78  );
79  const scalar maxIntSkew
80  (
81  readScalar(dict.lookup("maxInternalSkewness", true))
82  );
83  const scalar maxBounSkew
84  (
85  readScalar(dict.lookup("maxBoundarySkewness", true))
86  );
87  const scalar minWeight
88  (
89  readScalar(dict.lookup("minFaceWeight", true))
90  );
91  const scalar minVolRatio
92  (
93  readScalar(dict.lookup("minVolRatio", true))
94  );
95  const scalar minTwist
96  (
97  readScalar(dict.lookup("minTwist", true))
98  );
99  const scalar minTriangleTwist
100  (
101  readScalar(dict.lookup("minTriangleTwist", true))
102  );
103  const scalar minDet
104  (
105  readScalar(dict.lookup("minDeterminant", true))
106  );
107 
108  label nWrongFaces = 0;
109 
110  Info<< "Checking faces in error :" << endl;
111  //Pout.setf(ios_base::left);
112 
113  if (maxNonOrtho < 180.0-SMALL)
114  {
115  polyMeshGeometry::checkFaceDotProduct
116  (
117  report,
118  maxNonOrtho,
119  mesh,
120  mesh.cellCentres(),
121  mesh.faceAreas(),
122  checkFaces,
123  baffles,
124  &wrongFaces
125  );
126 
127  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
128 
129  Info<< " non-orthogonality > "
130  << setw(3) << maxNonOrtho
131  << " degrees : "
132  << nNewWrongFaces-nWrongFaces << endl;
133 
134  nWrongFaces = nNewWrongFaces;
135  }
136 
137  if (minVol > -GREAT)
138  {
139  polyMeshGeometry::checkFacePyramids
140  (
141  report,
142  minVol,
143  mesh,
144  mesh.cellCentres(),
145  mesh.points(),
146  checkFaces,
147  baffles,
148  &wrongFaces
149  );
150 
151  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
152 
153  Info<< " faces with face pyramid volume < "
154  << setw(5) << minVol << " : "
155  << nNewWrongFaces-nWrongFaces << endl;
156 
157  nWrongFaces = nNewWrongFaces;
158  }
159 
160  if (maxConcave < 180.0-SMALL)
161  {
162  polyMeshGeometry::checkFaceAngles
163  (
164  report,
165  maxConcave,
166  mesh,
167  mesh.faceAreas(),
168  mesh.points(),
169  checkFaces,
170  &wrongFaces
171  );
172 
173  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
174 
175  Info<< " faces with concavity > "
176  << setw(3) << maxConcave
177  << " degrees : "
178  << nNewWrongFaces-nWrongFaces << endl;
179 
180  nWrongFaces = nNewWrongFaces;
181  }
182 
183  if (minArea > -SMALL)
184  {
185  polyMeshGeometry::checkFaceArea
186  (
187  report,
188  minArea,
189  mesh,
190  mesh.faceAreas(),
191  checkFaces,
192  &wrongFaces
193  );
194 
195  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
196 
197  Info<< " faces with area < "
198  << setw(5) << minArea
199  << " m^2 : "
200  << nNewWrongFaces-nWrongFaces << endl;
201 
202  nWrongFaces = nNewWrongFaces;
203  }
204 
205  if (maxIntSkew > 0 || maxBounSkew > 0)
206  {
207  polyMeshGeometry::checkFaceSkewness
208  (
209  report,
210  maxIntSkew,
211  maxBounSkew,
212  mesh,
213  mesh.cellCentres(),
214  mesh.faceCentres(),
215  mesh.faceAreas(),
216  checkFaces,
217  baffles,
218  &wrongFaces
219  );
220 
221  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
222 
223  Info<< " faces with skewness > "
224  << setw(3) << maxIntSkew
225  << " (internal) or " << setw(3) << maxBounSkew
226  << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
227 
228  nWrongFaces = nNewWrongFaces;
229  }
230 
231  if (minWeight >= 0 && minWeight < 1)
232  {
233  polyMeshGeometry::checkFaceWeights
234  (
235  report,
236  minWeight,
237  mesh,
238  mesh.cellCentres(),
239  mesh.faceCentres(),
240  mesh.faceAreas(),
241  checkFaces,
242  baffles,
243  &wrongFaces
244  );
245 
246  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
247 
248  Info<< " faces with interpolation weights (0..1) < "
249  << setw(5) << minWeight
250  << " : "
251  << nNewWrongFaces-nWrongFaces << endl;
252 
253  nWrongFaces = nNewWrongFaces;
254  }
255 
256  if (minVolRatio >= 0)
257  {
258  polyMeshGeometry::checkVolRatio
259  (
260  report,
261  minVolRatio,
262  mesh,
263  mesh.cellVolumes(),
264  checkFaces,
265  baffles,
266  &wrongFaces
267  );
268 
269  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
270 
271  Info<< " faces with volume ratio of neighbour cells < "
272  << setw(5) << minVolRatio
273  << " : "
274  << nNewWrongFaces-nWrongFaces << endl;
275 
276  nWrongFaces = nNewWrongFaces;
277  }
278 
279  if (minTwist > -1)
280  {
281  //Pout<< "Checking face twist: dot product of face normal "
282  // << "with face triangle normals" << endl;
283  polyMeshGeometry::checkFaceTwist
284  (
285  report,
286  minTwist,
287  mesh,
288  mesh.cellCentres(),
289  mesh.faceAreas(),
290  mesh.faceCentres(),
291  mesh.points(),
292  checkFaces,
293  &wrongFaces
294  );
295 
296  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
297 
298  Info<< " faces with face twist < "
299  << setw(5) << minTwist
300  << " : "
301  << nNewWrongFaces-nWrongFaces << endl;
302 
303  nWrongFaces = nNewWrongFaces;
304  }
305 
306  if (minTriangleTwist > -1)
307  {
308  //Pout<< "Checking triangle twist: dot product of consecutive triangle"
309  // << " normals resulting from face-centre decomposition" << endl;
310  polyMeshGeometry::checkTriangleTwist
311  (
312  report,
313  minTriangleTwist,
314  mesh,
315  mesh.faceAreas(),
316  mesh.faceCentres(),
317  mesh.points(),
318  checkFaces,
319  &wrongFaces
320  );
321 
322  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
323 
324  Info<< " faces with triangle twist < "
325  << setw(5) << minTriangleTwist
326  << " : "
327  << nNewWrongFaces-nWrongFaces << endl;
328 
329  nWrongFaces = nNewWrongFaces;
330  }
331 
332  if (minDet > -1)
333  {
334  polyMeshGeometry::checkCellDeterminant
335  (
336  report,
337  minDet,
338  mesh,
339  mesh.faceAreas(),
340  checkFaces,
341  polyMeshGeometry::affectedCells(mesh, checkFaces),
342  &wrongFaces
343  );
344 
345  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
346 
347  Info<< " faces on cells with determinant < "
348  << setw(5) << minDet << " : "
349  << nNewWrongFaces-nWrongFaces << endl;
350 
351  nWrongFaces = nNewWrongFaces;
352  }
353 
354  //Pout.setf(ios_base::right);
355 
356  return nWrongFaces > 0;
357 }
358 
359 
361 (
362  const bool report,
363  const polyMesh& mesh,
364  const dictionary& dict,
365  labelHashSet& wrongFaces
366 )
367 {
368  return checkMesh
369  (
370  report,
371  mesh,
372  dict,
373  identity(mesh.nFaces()),
374  wrongFaces
375  );
376 }
377 
379 (
380  const bool report,
381  const dictionary& dict,
382  const polyMeshGeometry& meshGeom,
383  const labelList& checkFaces,
384  labelHashSet& wrongFaces
385 )
386 {
387  List<labelPair> emptyBaffles;
388 
389  return checkMesh
390  (
391  report,
392  dict,
393  meshGeom,
394  checkFaces,
395  emptyBaffles,
396  wrongFaces
397  );
398 }
399 
400 
402 (
403  const bool report,
404  const dictionary& dict,
405  const polyMeshGeometry& meshGeom,
406  const labelList& checkFaces,
407  const List<labelPair>& baffles,
408  labelHashSet& wrongFaces
409 )
410 {
411  const scalar maxNonOrtho
412  (
413  readScalar(dict.lookup("maxNonOrtho", true))
414  );
415  const scalar minVol
416  (
417  readScalar(dict.lookup("minVol", true))
418  );
419  const scalar maxConcave
420  (
421  readScalar(dict.lookup("maxConcave", true))
422  );
423  const scalar minArea
424  (
425  readScalar(dict.lookup("minArea", true))
426  );
427  const scalar maxIntSkew
428  (
429  readScalar(dict.lookup("maxInternalSkewness", true))
430  );
431  const scalar maxBounSkew
432  (
433  readScalar(dict.lookup("maxBoundarySkewness", true))
434  );
435  const scalar minWeight
436  (
437  readScalar(dict.lookup("minFaceWeight", true))
438  );
439  const scalar minVolRatio
440  (
441  readScalar(dict.lookup("minVolRatio", true))
442  );
443  const scalar minTwist
444  (
445  readScalar(dict.lookup("minTwist", true))
446  );
447  const scalar minTriangleTwist
448  (
449  readScalar(dict.lookup("minTriangleTwist", true))
450  );
451  const scalar minDet
452  (
453  readScalar(dict.lookup("minDeterminant", true))
454  );
455 
456  label nWrongFaces = 0;
457 
458  Info<< "Checking faces in error :" << endl;
459  //Pout.setf(ios_base::left);
460 
461  if (maxNonOrtho < 180.0-SMALL)
462  {
463  meshGeom.checkFaceDotProduct
464  (
465  report,
466  maxNonOrtho,
467  checkFaces,
468  baffles,
469  &wrongFaces
470  );
471 
472  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
473 
474  Info<< " non-orthogonality > "
475  << setw(3) << maxNonOrtho
476  << " degrees : "
477  << nNewWrongFaces-nWrongFaces << endl;
478 
479  nWrongFaces = nNewWrongFaces;
480  }
481 
482  if (minVol > -GREAT)
483  {
484  meshGeom.checkFacePyramids
485  (
486  report,
487  minVol,
488  meshGeom.mesh().points(),
489  checkFaces,
490  baffles,
491  &wrongFaces
492  );
493 
494  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
495 
496  Info<< " faces with face pyramid volume < "
497  << setw(5) << minVol << " : "
498  << nNewWrongFaces-nWrongFaces << endl;
499 
500  nWrongFaces = nNewWrongFaces;
501  }
502 
503  if (maxConcave < 180.0-SMALL)
504  {
505  meshGeom.checkFaceAngles
506  (
507  report,
508  maxConcave,
509  meshGeom.mesh().points(),
510  checkFaces,
511  &wrongFaces
512  );
513 
514  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
515 
516  Info<< " faces with concavity > "
517  << setw(3) << maxConcave
518  << " degrees : "
519  << nNewWrongFaces-nWrongFaces << endl;
520 
521  nWrongFaces = nNewWrongFaces;
522  }
523 
524  if (minArea > -SMALL)
525  {
526  meshGeom.checkFaceArea(report, minArea, checkFaces, &wrongFaces);
527 
528  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
529 
530  Info<< " faces with area < "
531  << setw(5) << minArea
532  << " m^2 : "
533  << nNewWrongFaces-nWrongFaces << endl;
534 
535  nWrongFaces = nNewWrongFaces;
536  }
537 
538  if (maxIntSkew > 0 || maxBounSkew > 0)
539  {
540  meshGeom.checkFaceSkewness
541  (
542  report,
543  maxIntSkew,
544  maxBounSkew,
545  checkFaces,
546  baffles,
547  &wrongFaces
548  );
549 
550  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
551 
552  Info<< " faces with skewness > "
553  << setw(3) << maxIntSkew
554  << " (internal) or " << setw(3) << maxBounSkew
555  << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
556 
557  nWrongFaces = nNewWrongFaces;
558  }
559 
560  if (minWeight >= 0 && minWeight < 1)
561  {
562  meshGeom.checkFaceWeights
563  (
564  report,
565  minWeight,
566  checkFaces,
567  baffles,
568  &wrongFaces
569  );
570 
571  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
572 
573  Info<< " faces with interpolation weights (0..1) < "
574  << setw(5) << minWeight
575  << " : "
576  << nNewWrongFaces-nWrongFaces << endl;
577 
578  nWrongFaces = nNewWrongFaces;
579  }
580 
581  if (minVolRatio >= 0)
582  {
583  meshGeom.checkVolRatio
584  (
585  report,
586  minVolRatio,
587  checkFaces,
588  baffles,
589  &wrongFaces
590  );
591 
592  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
593 
594  Info<< " faces with volume ratio of neighbour cells < "
595  << setw(5) << minVolRatio
596  << " : "
597  << nNewWrongFaces-nWrongFaces << endl;
598 
599  nWrongFaces = nNewWrongFaces;
600  }
601 
602  if (minTwist > -1)
603  {
604  //Pout<< "Checking face twist: dot product of face normal "
605  // << "with face triangle normals" << endl;
606  meshGeom.checkFaceTwist
607  (
608  report,
609  minTwist,
610  meshGeom.mesh().points(),
611  checkFaces,
612  &wrongFaces
613  );
614 
615  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
616 
617  Info<< " faces with face twist < "
618  << setw(5) << minTwist
619  << " : "
620  << nNewWrongFaces-nWrongFaces << endl;
621 
622  nWrongFaces = nNewWrongFaces;
623  }
624 
625  if (minTriangleTwist > -1)
626  {
627  //Pout<< "Checking triangle twist: dot product of consecutive triangle"
628  // << " normals resulting from face-centre decomposition" << endl;
629  meshGeom.checkTriangleTwist
630  (
631  report,
632  minTriangleTwist,
633  meshGeom.mesh().points(),
634  checkFaces,
635  &wrongFaces
636  );
637 
638  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
639 
640  Info<< " faces with triangle twist < "
641  << setw(5) << minTriangleTwist
642  << " : "
643  << nNewWrongFaces-nWrongFaces << endl;
644 
645  nWrongFaces = nNewWrongFaces;
646  }
647 
648  if (minDet > -1)
649  {
650  meshGeom.checkCellDeterminant
651  (
652  report,
653  minDet,
654  checkFaces,
655  meshGeom.affectedCells(meshGeom.mesh(), checkFaces),
656  &wrongFaces
657  );
658 
659  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
660 
661  Info<< " faces on cells with determinant < "
662  << setw(5) << minDet << " : "
663  << nNewWrongFaces-nWrongFaces << endl;
664 
665  nWrongFaces = nNewWrongFaces;
666  }
667 
668  //Pout.setf(ios_base::right);
669 
670  return nWrongFaces > 0;
671 }
672 
673 
674 // ************************ vim: set sw=4 sts=4 et: ************************ //