checkTopology.C 10.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
#include "checkTopology.H"
#include "polyMesh.H"
#include "Time.H"
#include "regionSplit.H"
#include "cellSet.H"
#include "faceSet.H"
#include "pointSet.H"
#include "IOmanip.H"

mattijs's avatar
mattijs committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
bool Foam::checkSync(const wordList& names)
{
    List<wordList> allNames(Pstream::nProcs());
    allNames[Pstream::myProcNo()] = names;
    Pstream::gatherList(allNames);

    bool hasError = false;

    for (label procI = 1; procI < allNames.size(); procI++)
    {
        if (allNames[procI] != allNames[0])
        {
            hasError = true;

            Info<< " ***Inconsistent zones across processors, "
                   "processor 0 has zones:" << allNames[0]
                << ", processor " << procI << " has zones:"
                << allNames[procI]
                << endl;
        }
    }
    return hasError;
}


35
36
37
38
39
40
Foam::label Foam::checkTopology
(
    const polyMesh& mesh,
    const bool allTopology,
    const bool allGeometry
)
41
42
43
{
    label noFailedChecks = 0;

mattijs's avatar
mattijs committed
44
    Info<< "Checking topology..." << endl;
45
46
47
48
49
50
51

    // Check if the boundary definition is unique
    mesh.boundaryMesh().checkDefinition(true);

    // Check if the boundary processor patches are correct
    mesh.boundaryMesh().checkParallelSync(true);

mattijs's avatar
mattijs committed
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
    // Check names of zones are equal
    if (checkSync(mesh.cellZones().names()))
    {
        noFailedChecks++;
    }
    if (checkSync(mesh.faceZones().names()))
    {
        noFailedChecks++;
    }
    if (checkSync(mesh.pointZones().names()))
    {
        noFailedChecks++;
    }

    // Check contents of faceZones consistent
    {
        forAll(mesh.faceZones(), zoneI)
        {
mattijs's avatar
mattijs committed
70
            if (mesh.faceZones()[zoneI].checkParallelSync(false))
mattijs's avatar
mattijs committed
71
            {
mattijs's avatar
mattijs committed
72
73
74
75
76
                Info<< " ***FaceZone " << mesh.faceZones()[zoneI].name()
                    << " is not correctly synchronised"
                    << " acrosss coupled boundaries."
                    << " (coupled faces both"
                    << " present in set but with opposite flipmap)" << endl;
mattijs's avatar
mattijs committed
77
78
79
80
81
                noFailedChecks++;
            }
        }
    }

82
83
84
85
86
87
    {
        pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
        if (mesh.checkPoints(true, &points))
        {
            noFailedChecks++;

mattijs's avatar
mattijs committed
88
89
90
            label nPoints = returnReduce(points.size(), sumOp<label>());

            Info<< "  <<Writing " << nPoints
91
                << " unused points to set " << points.name() << endl;
92
            points.instance() = mesh.pointsInstance();
93
94
95
96
97
98
99
100
101
            points.write();
        }
    }

    {
        faceSet faces(mesh, "upperTriangularFace", mesh.nFaces()/100);
        if (mesh.checkUpperTriangular(true, &faces))
        {
            noFailedChecks++;
102
        }
mattijs's avatar
mattijs committed
103
104
105
106

        label nFaces = returnReduce(faces.size(), sumOp<label>());

        if (nFaces > 0)
107
        {
mattijs's avatar
mattijs committed
108
            Info<< "  <<Writing " << nFaces
109
                << " unordered faces to set " << faces.name() << endl;
110
            faces.instance() = mesh.pointsInstance();
111
112
113
114
115
116
117
118
119
120
121
            faces.write();
        }
    }

    if (allTopology)
    {
        cellSet cells(mesh, "zipUpCells", mesh.nCells()/100);
        if (mesh.checkCellsZipUp(true, &cells))
        {
            noFailedChecks++;

mattijs's avatar
mattijs committed
122
123
124
            label nCells = returnReduce(cells.size(), sumOp<label>());

            Info<< "  <<Writing " << nCells
125
126
                << " cells with over used edges to set " << cells.name()
                << endl;
127
            cells.instance() = mesh.pointsInstance();
128
129
130
131
132
133
134
135
136
137
            cells.write();
        }
    }

    {
        faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
        if (mesh.checkFaceVertices(true, &faces))
        {
            noFailedChecks++;

mattijs's avatar
mattijs committed
138
139
140
            label nFaces = returnReduce(faces.size(), sumOp<label>());

            Info<< "  <<Writing " << nFaces
141
142
                << " faces with out-of-range or duplicate vertices to set "
                << faces.name() << endl;
143
            faces.instance() = mesh.pointsInstance();
144
145
146
147
148
149
150
151
152
153
154
            faces.write();
        }
    }

    if (allTopology)
    {
        faceSet faces(mesh, "edgeFaces", mesh.nFaces()/100);
        if (mesh.checkFaceFaces(true, &faces))
        {
            noFailedChecks++;

mattijs's avatar
mattijs committed
155
156
157
            label nFaces = returnReduce(faces.size(), sumOp<label>());

            Info<< "  <<Writing " << nFaces
158
159
                << " faces with incorrect edges to set " << faces.name()
                << endl;
160
            faces.instance() = mesh.pointsInstance();
161
162
163
164
            faces.write();
        }
    }

mattijs's avatar
mattijs committed
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
    if (allTopology)
    {
        labelList nInternalFaces(mesh.nCells(), 0);

        for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
        {
            nInternalFaces[mesh.faceOwner()[faceI]]++;
            nInternalFaces[mesh.faceNeighbour()[faceI]]++;
        }
        const polyBoundaryMesh& patches = mesh.boundaryMesh();
        forAll(patches, patchI)
        {
            if (patches[patchI].coupled())
            {
                const unallocLabelList& owners = patches[patchI].faceCells();

                forAll(owners, i)
                {
                    nInternalFaces[owners[i]]++;
                }
            }
        }

        faceSet oneCells(mesh, "oneInternalFaceCells", mesh.nCells()/100);
        faceSet twoCells(mesh, "twoInternalFacesCells", mesh.nCells()/100);

        forAll(nInternalFaces, cellI)
        {
            if (nInternalFaces[cellI] <= 1)
            {
                oneCells.insert(cellI);
            }
            else if (nInternalFaces[cellI] == 2)
            {
                twoCells.insert(cellI);
            }
        }

        label nOneCells = returnReduce(oneCells.size(), sumOp<label>());

        if (nOneCells > 0)
        {
            Info<< "  <<Writing " << nOneCells
                << " cells with with single non-boundary face to set "
                << oneCells.name()
                << endl;
211
            oneCells.instance() = mesh.pointsInstance();
mattijs's avatar
mattijs committed
212
213
214
215
216
217
218
219
220
221
222
            oneCells.write();
        }

        label nTwoCells = returnReduce(twoCells.size(), sumOp<label>());

        if (nTwoCells > 0)
        {
            Info<< "  <<Writing " << nTwoCells
                << " cells with with single non-boundary face to set "
                << twoCells.name()
                << endl;
223
            twoCells.instance() = mesh.pointsInstance();
mattijs's avatar
mattijs committed
224
225
226
227
            twoCells.write();
        }
    }

228
229
230
231
232
233
234
    {
        regionSplit rs(mesh);

        if (rs.nRegions() == 1)
        {
            Info<< "    Number of regions: " << rs.nRegions() << " (OK)."
                << endl;
235

236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
        }
        else
        {
            Info<< "   *Number of regions: " << rs.nRegions() << endl;

            Info<< "    The mesh has multiple regions which are not connected "
                   "by any face." << endl
                << "  <<Writing region information to "
                << mesh.time().timeName()/"cellToRegion"
                << endl;

            labelIOList ctr
            (
                IOobject
                (
                    "cellToRegion",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                rs
            );
            ctr.write();
        }
    }

    if (!Pstream::parRun())
    {
        Pout<< "\nChecking patch topology for multiply connected surfaces ..."
            << endl;

        const polyBoundaryMesh& patches = mesh.boundaryMesh();

        // Non-manifold points
        pointSet points
        (
            mesh,
            "nonManifoldPoints",
            mesh.nPoints()/100
        );

        Pout.setf(ios_base::left);

        Pout<< "    "
            << setw(20) << "Patch"
            << setw(9) << "Faces"
            << setw(9) << "Points"
284
285
286
287
288
289
            << setw(34) << "Surface topology";
        if (allGeometry)
        {
            Pout<< " Bounding box";
        }
        Pout<< endl;
290
291
292
293
294

        forAll(patches, patchI)
        {
            const polyPatch& pp = patches[patchI];

295
296
297
298
299
300
                Pout<< "    "
                    << setw(20) << pp.name()
                    << setw(9) << pp.size()
                    << setw(9) << pp.nPoints();


301
302
            primitivePatch::surfaceTopo pTyp = pp.surfaceType();

303
            if (pp.empty())
304
            {
305
                Pout<< setw(34) << "ok (empty)";
306
307
308
309
310
            }
            else if (pTyp == primitivePatch::MANIFOLD)
            {
                if (pp.checkPointManifold(true, &points))
                {
311
                    Pout<< setw(34) << "multiply connected (shared point)";
312
313
314
                }
                else
                {
315
                    Pout<< setw(34) << "ok (closed singly connected)";
316
317
318
319
320
321
322
323
324
325
326
                }

                // Add points on non-manifold edges to make set complete
                pp.checkTopology(false, &points);
            }
            else
            {
                pp.checkTopology(false, &points);

                if (pTyp == primitivePatch::OPEN)
                {
327
                    Pout<< setw(34) << "ok (non-closed singly connected)";
328
329
330
                }
                else
                {
331
332
333
334
335
336
337
338
339
                    Pout<< setw(34) << "multiply connected (shared edge)";
                }
            }

            if (allGeometry)
            {
                const pointField& pts = pp.points();
                const labelList& mp = pp.meshPoints();

340
                boundBox bb;   // zero-sized
341
342
343
344
345
346
347
348
349
350
351
                if (returnReduce(mp.size(), sumOp<label>()) > 0)
                {
                    bb.min() = pts[mp[0]];
                    bb.max() = pts[mp[0]];
                    for (label i = 1; i < mp.size(); i++)
                    {
                        bb.min() = min(bb.min(), pts[mp[i]]);
                        bb.max() = max(bb.max(), pts[mp[i]]);
                    }
                    reduce(bb.min(), minOp<vector>());
                    reduce(bb.max(), maxOp<vector>());
352
                }
353
                Pout<< ' ' << bb;
354
            }
355
            Pout<< endl;
356
357
        }

358
        if (points.size())
359
360
361
362
363
        {
            Pout<< "  <<Writing " << points.size()
                << " conflicting points to set "
                << points.name() << endl;

364
            points.instance() = mesh.pointsInstance();
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
            points.write();
        }

        //Pout.setf(ios_base::right);
    }

    // Force creation of all addressing if requested.
    // Errors will be reported as required
    if (allTopology)
    {
        mesh.cells();
        mesh.faces();
        mesh.edges();
        mesh.points();
        mesh.faceOwner();
        mesh.faceNeighbour();
        mesh.cellCells();
        mesh.edgeCells();
        mesh.pointCells();
        mesh.edgeFaces();
        mesh.pointFaces();
        mesh.cellEdges();
        mesh.faceEdges();
        mesh.pointEdges();
    }

    return noFailedChecks;
}