motionSmootherAlgoCheck.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-2020 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 "motionSmootherAlgo.H"
27 #include "polyMeshCheck.H"
28 #include "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  dict.lookup<scalar>("maxNonOrtho", true)
66  );
67  const scalar minVol
68  (
69  dict.lookup<scalar>("minVol", true)
70  );
71  const scalar minTetQuality
72  (
73  dict.lookup<scalar>("minTetQuality", true)
74  );
75  const scalar maxConcave
76  (
77  dict.lookup<scalar>("maxConcave", true)
78  );
79  const scalar minArea
80  (
81  dict.lookup<scalar>("minArea", true)
82  );
83  const scalar maxIntSkew
84  (
85  dict.lookup<scalar>("maxInternalSkewness", true)
86  );
87  const scalar maxBounSkew
88  (
89  dict.lookup<scalar>("maxBoundarySkewness", true)
90  );
91  const scalar minWeight
92  (
93  dict.lookup<scalar>("minFaceWeight", true)
94  );
95  const scalar minVolRatio
96  (
97  dict.lookup<scalar>("minVolRatio", true)
98  );
99  const scalar minTwist
100  (
101  dict.lookup<scalar>("minTwist", true)
102  );
103  const scalar minTriangleTwist
104  (
105  dict.lookup<scalar>("minTriangleTwist", true)
106  );
107  scalar minFaceFlatness = -1.0;
108  dict.readIfPresent("minFaceFlatness", minFaceFlatness, true);
109  const scalar minDet
110  (
111  dict.lookup<scalar>("minDeterminant", true)
112  );
113  label nWrongFaces = 0;
114 
115  Info<< "Checking faces in error :" << endl;
116  // Pout.setf(ios_base::left);
117 
118  if (maxNonOrtho < 180.0-small)
119  {
121  (
122  report,
123  maxNonOrtho,
124  mesh,
125  mesh.cellCentres(),
126  mesh.faceAreas(),
127  checkFaces,
128  baffles,
129  &wrongFaces
130  );
131 
132  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
133 
134  Info<< " non-orthogonality > "
135  << setw(3) << maxNonOrtho
136  << " degrees : "
137  << nNewWrongFaces-nWrongFaces << endl;
138 
139  nWrongFaces = nNewWrongFaces;
140  }
141 
142  if (minVol > -great)
143  {
145  (
146  report,
147  minVol,
148  mesh,
149  mesh.cellCentres(),
150  mesh.points(),
151  checkFaces,
152  baffles,
153  &wrongFaces
154  );
155 
156  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
157 
158  Info<< " faces with face pyramid volume < "
159  << setw(5) << minVol << " : "
160  << nNewWrongFaces-nWrongFaces << endl;
161 
162  nWrongFaces = nNewWrongFaces;
163  }
164 
165  if (minTetQuality > -great)
166  {
168  (
169  report,
170  minTetQuality,
171  mesh,
172  mesh.cellCentres(),
173  mesh.faceCentres(),
174  mesh.points(),
175  checkFaces,
176  baffles,
177  &wrongFaces
178  );
179 
180  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
181 
182  Info<< " faces with face-decomposition tet quality < "
183  << setw(5) << minTetQuality << " : "
184  << nNewWrongFaces-nWrongFaces << endl;
185 
186  nWrongFaces = nNewWrongFaces;
187  }
188 
189  if (maxConcave < 180.0-small)
190  {
192  (
193  report,
194  maxConcave,
195  mesh,
196  mesh.faceAreas(),
197  mesh.points(),
198  checkFaces,
199  &wrongFaces
200  );
201 
202  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
203 
204  Info<< " faces with concavity > "
205  << setw(3) << maxConcave
206  << " degrees : "
207  << nNewWrongFaces-nWrongFaces << endl;
208 
209  nWrongFaces = nNewWrongFaces;
210  }
211 
212  if (minArea > -small)
213  {
215  (
216  report,
217  minArea,
218  mesh,
219  mesh.faceAreas(),
220  checkFaces,
221  &wrongFaces
222  );
223 
224  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
225 
226  Info<< " faces with area < "
227  << setw(5) << minArea
228  << " m^2 : "
229  << nNewWrongFaces-nWrongFaces << endl;
230 
231  nWrongFaces = nNewWrongFaces;
232  }
233 
234  if (maxIntSkew > 0 || maxBounSkew > 0)
235  {
237  (
238  report,
239  maxIntSkew,
240  maxBounSkew,
241  mesh,
242  mesh.points(),
243  mesh.cellCentres(),
244  mesh.faceCentres(),
245  mesh.faceAreas(),
246  checkFaces,
247  baffles,
248  &wrongFaces
249  );
250 
251  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
252 
253  Info<< " faces with skewness > "
254  << setw(3) << maxIntSkew
255  << " (internal) or " << setw(3) << maxBounSkew
256  << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
257 
258  nWrongFaces = nNewWrongFaces;
259  }
260 
261  if (minWeight >= 0 && minWeight < 1)
262  {
264  (
265  report,
266  minWeight,
267  mesh,
268  mesh.cellCentres(),
269  mesh.faceCentres(),
270  mesh.faceAreas(),
271  checkFaces,
272  baffles,
273  &wrongFaces
274  );
275 
276  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
277 
278  Info<< " faces with interpolation weights (0..1) < "
279  << setw(5) << minWeight
280  << " : "
281  << nNewWrongFaces-nWrongFaces << endl;
282 
283  nWrongFaces = nNewWrongFaces;
284  }
285 
286  if (minVolRatio >= 0)
287  {
289  (
290  report,
291  minVolRatio,
292  mesh,
293  mesh.cellVolumes(),
294  checkFaces,
295  baffles,
296  &wrongFaces
297  );
298 
299  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
300 
301  Info<< " faces with volume ratio of neighbour cells < "
302  << setw(5) << minVolRatio
303  << " : "
304  << nNewWrongFaces-nWrongFaces << endl;
305 
306  nWrongFaces = nNewWrongFaces;
307  }
308 
309  if (minTwist > -1)
310  {
311  // Pout<< "Checking face twist: dot product of face normal "
312  // << "with face triangle normals" << endl;
314  (
315  report,
316  minTwist,
317  mesh,
318  mesh.cellCentres(),
319  mesh.faceAreas(),
320  mesh.faceCentres(),
321  mesh.points(),
322  checkFaces,
323  &wrongFaces
324  );
325 
326  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
327 
328  Info<< " faces with face twist < "
329  << setw(5) << minTwist
330  << " : "
331  << nNewWrongFaces-nWrongFaces << endl;
332 
333  nWrongFaces = nNewWrongFaces;
334  }
335 
336  if (minTriangleTwist > -1)
337  {
338  // Pout<< "Checking triangle twist: dot product of consecutive triangle"
339  // << " normals resulting from face-centre decomposition" << endl;
341  (
342  report,
343  minTriangleTwist,
344  mesh,
345  mesh.faceAreas(),
346  mesh.faceCentres(),
347  mesh.points(),
348  checkFaces,
349  &wrongFaces
350  );
351 
352  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
353 
354  Info<< " faces with triangle twist < "
355  << setw(5) << minTriangleTwist
356  << " : "
357  << nNewWrongFaces-nWrongFaces << endl;
358 
359  nWrongFaces = nNewWrongFaces;
360  }
361 
362  if (minFaceFlatness > -small)
363  {
365  (
366  report,
367  minFaceFlatness,
368  mesh,
369  mesh.faceAreas(),
370  mesh.faceCentres(),
371  mesh.points(),
372  checkFaces,
373  &wrongFaces
374  );
375 
376  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
377 
378  Info<< " faces with flatness < "
379  << setw(5) << minFaceFlatness
380  << " : "
381  << nNewWrongFaces-nWrongFaces << endl;
382 
383  nWrongFaces = nNewWrongFaces;
384  }
385 
386  if (minDet > -1)
387  {
389  (
390  report,
391  minDet,
392  mesh,
393  mesh.faceAreas(),
394  checkFaces,
395  &wrongFaces
396  );
397 
398  label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
399 
400  Info<< " faces on cells with determinant < "
401  << setw(5) << minDet << " : "
402  << nNewWrongFaces-nWrongFaces << endl;
403 
404  nWrongFaces = nNewWrongFaces;
405  }
406 
407  // Pout.setf(ios_base::right);
408 
409  return nWrongFaces > 0;
410 }
411 
412 
414 (
415  const bool report,
416  const polyMesh& mesh,
417  const dictionary& dict,
418  labelHashSet& wrongFaces
419 )
420 {
421  return checkMesh
422  (
423  report,
424  mesh,
425  dict,
426  identity(mesh.nFaces()),
427  wrongFaces
428  );
429 }
430 
431 
432 // ************************************************************************* //
static bool checkMesh(const bool report, const polyMesh &mesh, const dictionary &dict, labelHashSet &wrongFaces)
Check mesh with mesh settings in dict. Collects incorrect faces.
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.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
bool checkFaceArea(const bool report, const scalar minArea, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
Check for small faces.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
label nFaces() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
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 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.
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
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.
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).
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
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)
const vectorField & cellCentres() const
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.
Istream and Ostream manipulators taking arguments.
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.
const vectorField & faceCentres() const
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.
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.
const vectorField & faceAreas() const
messageStream Info
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
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:74
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
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)
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844
const scalarField & cellVolumes() const