checkTopology.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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 "checkTopology.H"
27 #include "polyMesh.H"
28 #include "Time.H"
29 #include "regionSplit.H"
30 #include "cellSet.H"
31 #include "faceSet.H"
32 #include "pointSet.H"
33 #include "IOmanip.H"
34 #include "emptyPolyPatch.H"
35 #include "processorPolyPatch.H"
36 #include "surfaceWriter.H"
37 #include "checkTools.H"
38 
39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
40 
42 (
43  const polyMesh& mesh,
44  const bool allTopology,
45  const bool allGeometry,
46  const autoPtr<surfaceWriter>& writer
47 )
48 {
49  label noFailedChecks = 0;
50 
51  Info<< "Checking topology..." << endl;
52 
53  // Check if the boundary definition is unique
54  mesh.boundaryMesh().checkDefinition(true);
55 
56  // Check that empty patches cover all sides of the mesh
57  {
58  label nEmpty = 0;
59  forAll(mesh.boundaryMesh(), patchi)
60  {
61  if (isA<emptyPolyPatch>(mesh.boundaryMesh()[patchi]))
62  {
63  nEmpty += mesh.boundaryMesh()[patchi].size();
64  }
65  }
66  reduce(nEmpty, sumOp<label>());
67  label nTotCells = returnReduce(mesh.cells().size(), sumOp<label>());
68 
69  // These are actually warnings, not errors.
70  if (nTotCells && (nEmpty % nTotCells))
71  {
72  Info<< " ***Total number of faces on empty patches"
73  << " is not divisible by the number of cells in the mesh."
74  << " Hence this mesh is not 1D or 2D."
75  << endl;
76  }
77  }
78 
79  // Check if the boundary processor patches are correct
80  mesh.boundaryMesh().checkParallelSync(true);
81 
82  // Check names of zones are equal
83  mesh.cellZones().checkDefinition(true);
84  if (mesh.cellZones().checkParallelSync(true))
85  {
86  noFailedChecks++;
87  }
88  mesh.faceZones().checkDefinition(true);
89  if (mesh.faceZones().checkParallelSync(true))
90  {
91  noFailedChecks++;
92  }
93  mesh.pointZones().checkDefinition(true);
94  if (mesh.pointZones().checkParallelSync(true))
95  {
96  noFailedChecks++;
97  }
98 
99 
100  {
101  cellSet cells(mesh, "illegalCells", mesh.nCells()/100);
102 
103  forAll(mesh.cells(), celli)
104  {
105  const cell& cFaces = mesh.cells()[celli];
106 
107  if (cFaces.size() <= 3)
108  {
109  cells.insert(celli);
110  }
111  forAll(cFaces, i)
112  {
113  if (cFaces[i] < 0 || cFaces[i] >= mesh.nFaces())
114  {
115  cells.insert(celli);
116  break;
117  }
118  }
119  }
120  label nCells = returnReduce(cells.size(), sumOp<label>());
121 
122  if (nCells > 0)
123  {
124  Info<< " Illegal cells (less than 4 faces or out of range faces)"
125  << " found, number of cells: " << nCells << endl;
126  noFailedChecks++;
127 
128  Info<< " <<Writing " << nCells
129  << " illegal cells to set " << cells.name() << endl;
130  cells.instance() = mesh.pointsInstance();
131  cells.write();
132  if (writer.valid())
133  {
134  mergeAndWrite(writer(), cells);
135  }
136 
137  }
138  else
139  {
140  Info<< " Cell to face addressing OK." << endl;
141  }
142  }
143 
144 
145  {
146  pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
147  if (mesh.checkPoints(true, &points))
148  {
149  noFailedChecks++;
150 
151  label nPoints = returnReduce(points.size(), sumOp<label>());
152 
153  Info<< " <<Writing " << nPoints
154  << " unused points to set " << points.name() << endl;
155  points.instance() = mesh.pointsInstance();
156  points.write();
157  }
158  }
159 
160  {
161  faceSet faces(mesh, "upperTriangularFace", mesh.nFaces()/100);
162  if (mesh.checkUpperTriangular(true, &faces))
163  {
164  noFailedChecks++;
165  }
166 
167  label nFaces = returnReduce(faces.size(), sumOp<label>());
168 
169  if (nFaces > 0)
170  {
171  Info<< " <<Writing " << nFaces
172  << " unordered faces to set " << faces.name() << endl;
173  faces.instance() = mesh.pointsInstance();
174  faces.write();
175  if (writer.valid())
176  {
177  mergeAndWrite(writer(), faces);
178  }
179  }
180  }
181 
182  {
183  faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
184  if (mesh.checkFaceVertices(true, &faces))
185  {
186  noFailedChecks++;
187 
188  label nFaces = returnReduce(faces.size(), sumOp<label>());
189 
190  Info<< " <<Writing " << nFaces
191  << " faces with out-of-range or duplicate vertices to set "
192  << faces.name() << endl;
193  faces.instance() = mesh.pointsInstance();
194  faces.write();
195  if (writer.valid())
196  {
197  mergeAndWrite(writer(), faces);
198  }
199  }
200  }
201 
202  if (allTopology)
203  {
204  cellSet cells(mesh, "zipUpCells", mesh.nCells()/100);
205  if (mesh.checkCellsZipUp(true, &cells))
206  {
207  noFailedChecks++;
208 
209  label nCells = returnReduce(cells.size(), sumOp<label>());
210 
211  Info<< " <<Writing " << nCells
212  << " cells with over used edges to set " << cells.name()
213  << endl;
214  cells.instance() = mesh.pointsInstance();
215  cells.write();
216  if (writer.valid())
217  {
218  mergeAndWrite(writer(), cells);
219  }
220 
221  }
222  }
223 
224  if (allTopology)
225  {
226  faceSet faces(mesh, "edgeFaces", mesh.nFaces()/100);
227  if (mesh.checkFaceFaces(true, &faces))
228  {
229  noFailedChecks++;
230  }
231 
232  label nFaces = returnReduce(faces.size(), sumOp<label>());
233  if (nFaces > 0)
234  {
235  Info<< " <<Writing " << nFaces
236  << " faces with non-standard edge connectivity to set "
237  << faces.name() << endl;
238  faces.instance() = mesh.pointsInstance();
239  faces.write();
240  if (writer.valid())
241  {
242  mergeAndWrite(writer(), faces);
243  }
244  }
245  }
246 
247  if (allTopology)
248  {
249  labelList nInternalFaces(mesh.nCells(), 0);
250 
251  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
252  {
253  nInternalFaces[mesh.faceOwner()[facei]]++;
254  nInternalFaces[mesh.faceNeighbour()[facei]]++;
255  }
256  const polyBoundaryMesh& patches = mesh.boundaryMesh();
257  forAll(patches, patchi)
258  {
259  if (patches[patchi].coupled())
260  {
261  const labelUList& owners = patches[patchi].faceCells();
262 
263  forAll(owners, i)
264  {
265  nInternalFaces[owners[i]]++;
266  }
267  }
268  }
269 
270  cellSet oneCells(mesh, "oneInternalFaceCells", mesh.nCells()/100);
271  cellSet twoCells(mesh, "twoInternalFacesCells", mesh.nCells()/100);
272 
273  forAll(nInternalFaces, celli)
274  {
275  if (nInternalFaces[celli] <= 1)
276  {
277  oneCells.insert(celli);
278  }
279  else if (nInternalFaces[celli] == 2)
280  {
281  twoCells.insert(celli);
282  }
283  }
284 
285  label nOneCells = returnReduce(oneCells.size(), sumOp<label>());
286 
287  if (nOneCells > 0)
288  {
289  Info<< " <<Writing " << nOneCells
290  << " cells with zero or one non-boundary face to set "
291  << oneCells.name()
292  << endl;
293  oneCells.instance() = mesh.pointsInstance();
294  oneCells.write();
295  if (writer.valid())
296  {
297  mergeAndWrite(writer(), oneCells);
298  }
299  }
300 
301  label nTwoCells = returnReduce(twoCells.size(), sumOp<label>());
302 
303  if (nTwoCells > 0)
304  {
305  Info<< " <<Writing " << nTwoCells
306  << " cells with two non-boundary faces to set "
307  << twoCells.name()
308  << endl;
309  twoCells.instance() = mesh.pointsInstance();
310  twoCells.write();
311  if (writer.valid())
312  {
313  mergeAndWrite(writer(), twoCells);
314  }
315  }
316  }
317 
318  {
319  regionSplit rs(mesh);
320 
321  if (rs.nRegions() <= 1)
322  {
323  Info<< " Number of regions: " << rs.nRegions() << " (OK)."
324  << endl;
325 
326  }
327  else
328  {
329  Info<< " *Number of regions: " << rs.nRegions() << endl;
330 
331  Info<< " The mesh has multiple regions which are not connected "
332  "by any face." << endl
333  << " <<Writing region information to "
334  << mesh.time().timeName()/"cellToRegion"
335  << endl;
336 
337  labelIOList ctr
338  (
339  IOobject
340  (
341  "cellToRegion",
342  mesh.time().timeName(),
343  mesh,
344  IOobject::NO_READ,
345  IOobject::NO_WRITE
346  ),
347  rs
348  );
349  ctr.write();
350 
351 
352 
353  // Is region disconnected
354  boolList regionDisconnected(rs.nRegions(), true);
355  if (allTopology)
356  {
357  // -1 : not assigned
358  // -2 : multiple regions
359  // >= 0 : single region
360  labelList pointToRegion(mesh.nPoints(), -1);
361 
362  for
363  (
364  label faceI = mesh.nInternalFaces();
365  faceI < mesh.nFaces();
366  faceI++
367  )
368  {
369  label regionI = rs[mesh.faceOwner()[faceI]];
370  const face& f = mesh.faces()[faceI];
371  forAll(f, fp)
372  {
373  label& pRegion = pointToRegion[f[fp]];
374  if (pRegion == -1)
375  {
376  pRegion = regionI;
377  }
378  else if (pRegion == -2)
379  {
380  // Already marked
381  regionDisconnected[regionI] = false;
382  }
383  else if (pRegion != regionI)
384  {
385  // Multiple regions
386  regionDisconnected[regionI] = false;
387  regionDisconnected[pRegion] = false;
388  pRegion = -2;
389  }
390  }
391  }
392 
393  Pstream::listCombineGather(regionDisconnected, andEqOp<bool>());
394  Pstream::listCombineScatter(regionDisconnected);
395  }
396 
397 
398 
399  // write cellSet for each region
400  PtrList<cellSet> cellRegions(rs.nRegions());
401  for (label i = 0; i < rs.nRegions(); i++)
402  {
403  cellRegions.set
404  (
405  i,
406  new cellSet
407  (
408  mesh,
409  "region" + Foam::name(i),
410  mesh.nCells()/100
411  )
412  );
413  }
414 
415  forAll(rs, i)
416  {
417  cellRegions[rs[i]].insert(i);
418  }
419 
420  for (label i = 0; i < rs.nRegions(); i++)
421  {
422  Info<< " <<Writing region " << i;
423  if (allTopology)
424  {
425  if (regionDisconnected[i])
426  {
427  Info<< " (fully disconnected)";
428  }
429  else
430  {
431  Info<< " (point connected)";
432  }
433  }
434  Info<< " with "
435  << returnReduce(cellRegions[i].size(), sumOp<scalar>())
436  << " cells to cellSet " << cellRegions[i].name() << endl;
437 
438  cellRegions[i].write();
439  }
440  }
441  }
442 
443 
444  {
445  if (!Pstream::parRun())
446  {
447  Info<< "\nChecking patch topology for multiply connected"
448  << " surfaces..." << endl;
449  }
450  else
451  {
452  Info<< "\nChecking basic patch addressing..." << endl;
453  }
454 
455 
456  const polyBoundaryMesh& patches = mesh.boundaryMesh();
457 
458  // Non-manifold points
459  pointSet points
460  (
461  mesh,
462  "nonManifoldPoints",
463  mesh.nPoints()/1000
464  );
465 
466  Pout.setf(ios_base::left);
467 
468  Info<< " "
469  << setw(20) << "Patch"
470  << setw(9) << "Faces"
471  << setw(9) << "Points";
472  if (!Pstream::parRun())
473  {
474  Info<< setw(34) << "Surface topology";
475  }
476  if (allGeometry)
477  {
478  Info<< " Bounding box";
479  }
480  Info<< endl;
481 
482  forAll(patches, patchi)
483  {
484  const polyPatch& pp = patches[patchi];
485 
486  if (!isA<processorPolyPatch>(pp))
487  {
488  Info<< " "
489  << setw(20) << pp.name()
490  << setw(9) << returnReduce(pp.size(), sumOp<label>())
491  << setw(9) << returnReduce(pp.nPoints(), sumOp<label>());
492 
493  if (!Pstream::parRun())
494  {
495  primitivePatch::surfaceTopo pTyp = pp.surfaceType();
496 
497  if (pp.empty())
498  {
499  Info<< setw(34) << "ok (empty)";
500  }
501  else if (pTyp == primitivePatch::MANIFOLD)
502  {
503  if (pp.checkPointManifold(true, &points))
504  {
505  Info<< setw(34)
506  << "multiply connected (shared point)";
507  }
508  else
509  {
510  Info<< setw(34) << "ok (closed singly connected)";
511  }
512 
513  // Add points on non-manifold edges to make set complete
514  pp.checkTopology(false, &points);
515  }
516  else
517  {
518  pp.checkTopology(false, &points);
519 
520  if (pTyp == primitivePatch::OPEN)
521  {
522  Info<< setw(34)
523  << "ok (non-closed singly connected)";
524  }
525  else
526  {
527  Info<< setw(34)
528  << "multiply connected (shared edge)";
529  }
530  }
531  }
532 
533  if (allGeometry)
534  {
535  const pointField& pts = pp.points();
536  const labelList& mp = pp.meshPoints();
537 
538  if (returnReduce(mp.size(), sumOp<label>()) > 0)
539  {
540  boundBox bb(point::max, point::min);
541  forAll(mp, i)
542  {
543  bb.min() = min(bb.min(), pts[mp[i]]);
544  bb.max() = max(bb.max(), pts[mp[i]]);
545  }
546  reduce(bb.min(), minOp<vector>());
547  reduce(bb.max(), maxOp<vector>());
548  Info<< ' ' << bb;
549  }
550  }
551  Info<< endl;
552  }
553  }
554 
555  if (points.size())
556  {
557  Info<< " <<Writing " << returnReduce(points.size(), sumOp<label>())
558  << " conflicting points to set "
559  << points.name() << endl;
560 
561  points.instance() = mesh.pointsInstance();
562  points.write();
563  }
564 
565  //Info.setf(ios_base::right);
566  }
567 
568  // Force creation of all addressing if requested.
569  // Errors will be reported as required
570  if (allTopology)
571  {
572  mesh.cells();
573  mesh.faces();
574  mesh.edges();
575  mesh.points();
576  mesh.faceOwner();
577  mesh.faceNeighbour();
578  mesh.cellCells();
579  mesh.edgeCells();
580  mesh.pointCells();
581  mesh.edgeFaces();
582  mesh.pointFaces();
583  mesh.cellEdges();
584  mesh.faceEdges();
585  mesh.pointEdges();
586  }
587 
588  return noFailedChecks;
589 }
590 
591 
592 // ************************************************************************* //
label checkTopology(const polyMesh &, const bool, const bool, const autoPtr< surfaceWriter > &)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
ios_base::fmtflags setf(const ios_base::fmtflags f)
Set flags of stream.
Definition: IOstream.H:496
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
UList< label > labelUList
Definition: UList.H:64
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
dynamicFvMesh & mesh
const cellShapeList & cells
const pointField & points
List< label > labelList
A List of labels.
Definition: labelList.H:56
void mergeAndWrite(const polyMesh &mesh, const surfaceWriter &writer, const word &name, const indirectPrimitivePatch setPatch, const fileName &outputDir)
Generate merged surface on master and write. Needs input patch.
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
Istream and Ostream manipulators taking arguments.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
label patchi
messageStream Info
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
virtual Ostream & write(const token &)=0
Write next token to stream.
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
IOList< label > labelIOList
Label container classes.
Definition: labelIOList.H:42