applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

   1 #include "checkTopology.H"
   2 #include "polyMesh.H"
   3 #include "Time.H"
   4 #include "regionSplit.H"
   5 #include "cellSet.H"
   6 #include "faceSet.H"
   7 #include "pointSet.H"
   8 #include "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         faceSet oneCells(mesh, "oneInternalFaceCells", mesh.nCells()/100);
 185         faceSet 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                 boundBox bb;   // zero-sized
 335                 if (returnReduce(mp.size(), sumOp<label>()) > 0)
 336                 {
 337                     bb.min() = pts[mp[0]];
 338                     bb.max() = pts[mp[0]];
 339                     for (label i = 1; i < mp.size(); i++)
 340                     {
 341                         bb.min() = min(bb.min(), pts[mp[i]]);
 342                         bb.max() = max(bb.max(), pts[mp[i]]);
 343                     }
 344                     reduce(bb.min(), minOp<vector>());
 345                     reduce(bb.max(), maxOp<vector>());
 346                 }
 347                 Pout<< ' ' << bb;
 348             }
 349             Pout<< endl;
 350         }
 351
 352         if (points.size())
 353         {
 354             Pout<< "  <<Writing " << points.size()
 355                 << " conflicting points to set "
 356                 << points.name() << endl;
 357
 358             points.write();
 359         }
 360
 361         //Pout.setf(ios_base::right);
 362     }
 363
 364     // Force creation of all addressing if requested.
 365     // Errors will be reported as required
 366     if (allTopology)
 367     {
 368         mesh.cells();
 369         mesh.faces();
 370         mesh.edges();
 371         mesh.points();
 372         mesh.faceOwner();
 373         mesh.faceNeighbour();
 374         mesh.cellCells();
 375         mesh.edgeCells();
 376         mesh.pointCells();
 377         mesh.edgeFaces();
 378         mesh.pointFaces();
 379         mesh.cellEdges();
 380         mesh.faceEdges();
 381         mesh.pointEdges();
 382     }
 383
 384     return noFailedChecks;
 385 }