FreeFOAM The Cross-Platform CFD Toolkit
checkTopology.C
Go to the documentation of this file.
1 #include "checkTopology.H"
2 #include <OpenFOAM/polyMesh.H>
3 #include <OpenFOAM/Time.H>
5 #include <meshTools/cellSet.H>
6 #include <meshTools/faceSet.H>
7 #include <meshTools/pointSet.H>
8 #include <OpenFOAM/IOmanip.H>
9 
10 bool Foam::checkSync(const wordList& names)
11 {
12  List<wordList> allNames(Pstream::nProcs());
13  allNames[Pstream::myProcNo()] = names;
14  Pstream::gatherList(allNames);
15  Pstream::scatterList(allNames);
16 
17  bool hasError = false;
18 
19  for (label procI = 1; procI < allNames.size(); procI++)
20  {
21  if (allNames[procI] != allNames[0])
22  {
23  hasError = true;
24 
25  Info<< " ***Inconsistent zones across processors, "
26  "processor 0 has zones:" << allNames[0]
27  << ", processor " << procI << " has zones:"
28  << allNames[procI]
29  << endl;
30  }
31  }
32  return hasError;
33 }
34 
35 
36 Foam::label Foam::checkTopology
37 (
38  const polyMesh& mesh,
39  const bool allTopology,
40  const bool allGeometry
41 )
42 {
43  label noFailedChecks = 0;
44 
45  Info<< "Checking topology..." << endl;
46 
47  // Check if the boundary definition is unique
48  mesh.boundaryMesh().checkDefinition(true);
49 
50  // Check if the boundary processor patches are correct
51  mesh.boundaryMesh().checkParallelSync(true);
52 
53  // Check names of zones are equal
54  if (checkSync(mesh.cellZones().names()))
55  {
56  noFailedChecks++;
57  }
58  if (checkSync(mesh.faceZones().names()))
59  {
60  noFailedChecks++;
61  }
62  if (checkSync(mesh.pointZones().names()))
63  {
64  noFailedChecks++;
65  }
66 
67  // Check contents of faceZones consistent
68  {
69  forAll(mesh.faceZones(), zoneI)
70  {
71  if (mesh.faceZones()[zoneI].checkParallelSync(false))
72  {
73  Info<< " ***FaceZone " << mesh.faceZones()[zoneI].name()
74  << " is not correctly synchronised"
75  << " across coupled boundaries."
76  << " (coupled faces are either not both "
77  << " present in set or have same flipmap)" << endl;
78  noFailedChecks++;
79  }
80  }
81  }
82 
83  {
84  pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
85  if (mesh.checkPoints(true, &points))
86  {
87  noFailedChecks++;
88 
89  label nPoints = returnReduce(points.size(), sumOp<label>());
90 
91  Info<< " <<Writing " << nPoints
92  << " unused points to set " << points.name() << endl;
93  points.write();
94  }
95  }
96 
97  {
98  faceSet faces(mesh, "upperTriangularFace", mesh.nFaces()/100);
99  if (mesh.checkUpperTriangular(true, &faces))
100  {
101  noFailedChecks++;
102  }
103 
104  label nFaces = returnReduce(faces.size(), sumOp<label>());
105 
106  if (nFaces > 0)
107  {
108  Info<< " <<Writing " << nFaces
109  << " unordered faces to set " << faces.name() << endl;
110  faces.write();
111  }
112  }
113 
114  if (allTopology)
115  {
116  cellSet cells(mesh, "zipUpCells", mesh.nCells()/100);
117  if (mesh.checkCellsZipUp(true, &cells))
118  {
119  noFailedChecks++;
120 
121  label nCells = returnReduce(cells.size(), sumOp<label>());
122 
123  Info<< " <<Writing " << nCells
124  << " cells with over used edges to set " << cells.name()
125  << endl;
126  cells.write();
127  }
128  }
129 
130  {
131  faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
132  if (mesh.checkFaceVertices(true, &faces))
133  {
134  noFailedChecks++;
135 
136  label nFaces = returnReduce(faces.size(), sumOp<label>());
137 
138  Info<< " <<Writing " << nFaces
139  << " faces with out-of-range or duplicate vertices to set "
140  << faces.name() << endl;
141  faces.write();
142  }
143  }
144 
145  if (allTopology)
146  {
147  faceSet faces(mesh, "edgeFaces", mesh.nFaces()/100);
148  if (mesh.checkFaceFaces(true, &faces))
149  {
150  noFailedChecks++;
151 
152  label nFaces = returnReduce(faces.size(), sumOp<label>());
153 
154  Info<< " <<Writing " << nFaces
155  << " faces with incorrect edges to set " << faces.name()
156  << endl;
157  faces.write();
158  }
159  }
160 
161  if (allTopology)
162  {
163  labelList nInternalFaces(mesh.nCells(), 0);
164 
165  for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
166  {
167  nInternalFaces[mesh.faceOwner()[faceI]]++;
168  nInternalFaces[mesh.faceNeighbour()[faceI]]++;
169  }
170  const polyBoundaryMesh& patches = mesh.boundaryMesh();
171  forAll(patches, patchI)
172  {
173  if (patches[patchI].coupled())
174  {
175  const unallocLabelList& owners = patches[patchI].faceCells();
176 
177  forAll(owners, i)
178  {
179  nInternalFaces[owners[i]]++;
180  }
181  }
182  }
183 
184  cellSet oneCells(mesh, "oneInternalFaceCells", mesh.nCells()/100);
185  cellSet twoCells(mesh, "twoInternalFacesCells", mesh.nCells()/100);
186 
187  forAll(nInternalFaces, cellI)
188  {
189  if (nInternalFaces[cellI] <= 1)
190  {
191  oneCells.insert(cellI);
192  }
193  else if (nInternalFaces[cellI] == 2)
194  {
195  twoCells.insert(cellI);
196  }
197  }
198 
199  label nOneCells = returnReduce(oneCells.size(), sumOp<label>());
200 
201  if (nOneCells > 0)
202  {
203  Info<< " <<Writing " << nOneCells
204  << " cells with with single non-boundary face to set "
205  << oneCells.name()
206  << endl;
207  oneCells.write();
208  }
209 
210  label nTwoCells = returnReduce(twoCells.size(), sumOp<label>());
211 
212  if (nTwoCells > 0)
213  {
214  Info<< " <<Writing " << nTwoCells
215  << " cells with with single non-boundary face to set "
216  << twoCells.name()
217  << endl;
218  twoCells.write();
219  }
220  }
221 
222  {
223  regionSplit rs(mesh);
224 
225  if (rs.nRegions() == 1)
226  {
227  Info<< " Number of regions: " << rs.nRegions() << " (OK)."
228  << endl;
229 
230  }
231  else
232  {
233  Info<< " *Number of regions: " << rs.nRegions() << endl;
234 
235  Info<< " The mesh has multiple regions which are not connected "
236  "by any face." << endl
237  << " <<Writing region information to "
238  << mesh.time().timeName()/"cellToRegion"
239  << endl;
240 
241  labelIOList ctr
242  (
243  IOobject
244  (
245  "cellToRegion",
246  mesh.time().timeName(),
247  mesh,
248  IOobject::NO_READ,
249  IOobject::NO_WRITE
250  ),
251  rs
252  );
253  ctr.write();
254  }
255  }
256 
257  if (!Pstream::parRun())
258  {
259  Pout<< "\nChecking patch topology for multiply connected surfaces ..."
260  << endl;
261 
262  const polyBoundaryMesh& patches = mesh.boundaryMesh();
263 
264  // Non-manifold points
265  pointSet points
266  (
267  mesh,
268  "nonManifoldPoints",
269  mesh.nPoints()/100
270  );
271 
272  Pout.setf(ios_base::left);
273 
274  Pout<< " "
275  << setw(20) << "Patch"
276  << setw(9) << "Faces"
277  << setw(9) << "Points"
278  << setw(34) << "Surface topology";
279  if (allGeometry)
280  {
281  Pout<< " Bounding box";
282  }
283  Pout<< endl;
284 
285  forAll(patches, patchI)
286  {
287  const polyPatch& pp = patches[patchI];
288 
289  Pout<< " "
290  << setw(20) << pp.name()
291  << setw(9) << pp.size()
292  << setw(9) << pp.nPoints();
293 
294 
295  primitivePatch::surfaceTopo pTyp = pp.surfaceType();
296 
297  if (pp.empty())
298  {
299  Pout<< setw(34) << "ok (empty)";
300  }
301  else if (pTyp == primitivePatch::MANIFOLD)
302  {
303  if (pp.checkPointManifold(true, &points))
304  {
305  Pout<< setw(34) << "multiply connected (shared point)";
306  }
307  else
308  {
309  Pout<< setw(34) << "ok (closed singly connected)";
310  }
311 
312  // Add points on non-manifold edges to make set complete
313  pp.checkTopology(false, &points);
314  }
315  else
316  {
317  pp.checkTopology(false, &points);
318 
319  if (pTyp == primitivePatch::OPEN)
320  {
321  Pout<< setw(34) << "ok (non-closed singly connected)";
322  }
323  else
324  {
325  Pout<< setw(34) << "multiply connected (shared edge)";
326  }
327  }
328 
329  if (allGeometry)
330  {
331  const pointField& pts = pp.points();
332  const labelList& mp = pp.meshPoints();
333 
334  if (returnReduce(mp.size(), sumOp<label>()) > 0)
335  {
336  boundBox bb(point::max, point::min);
337  forAll(mp, i)
338  {
339  bb.min() = min(bb.min(), pts[mp[i]]);
340  bb.max() = max(bb.max(), pts[mp[i]]);
341  }
342  reduce(bb.min(), minOp<vector>());
343  reduce(bb.max(), maxOp<vector>());
344  Pout<< ' ' << bb;
345  }
346  }
347  Pout<< endl;
348  }
349 
350  if (points.size())
351  {
352  Pout<< " <<Writing " << points.size()
353  << " conflicting points to set "
354  << points.name() << endl;
355 
356  points.write();
357  }
358 
359  //Pout.setf(ios_base::right);
360  }
361 
362  // Force creation of all addressing if requested.
363  // Errors will be reported as required
364  if (allTopology)
365  {
366  mesh.cells();
367  mesh.faces();
368  mesh.edges();
369  mesh.points();
370  mesh.faceOwner();
371  mesh.faceNeighbour();
372  mesh.cellCells();
373  mesh.edgeCells();
374  mesh.pointCells();
375  mesh.edgeFaces();
376  mesh.pointFaces();
377  mesh.cellEdges();
378  mesh.faceEdges();
379  mesh.pointEdges();
380  }
381 
382  return noFailedChecks;
383 }