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
bool Foam::checkSync(const wordList& names)
{
    List<wordList> allNames(Pstream::nProcs());
    allNames[Pstream::myProcNo()] = names;
    Pstream::gatherList(allNames);
15
    Pstream::scatterList(allNames);
mattijs's avatar
mattijs committed
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

    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;
}


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

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

    // 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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
    // 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
71
            if (mesh.faceZones()[zoneI].checkParallelSync(false))
mattijs's avatar
mattijs committed
72
            {
mattijs's avatar
mattijs committed
73
74
                Info<< " ***FaceZone " << mesh.faceZones()[zoneI].name()
                    << " is not correctly synchronised"
mattijs's avatar
mattijs committed
75
                    << " across coupled boundaries."
76
77
                    << " (coupled faces are either not both "
                    << " present in set or have same flipmap)" << endl;
mattijs's avatar
mattijs committed
78
79
80
81
82
                noFailedChecks++;
            }
        }
    }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

mattijs's avatar
mattijs committed
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
211
    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;
212
            oneCells.instance() = mesh.pointsInstance();
mattijs's avatar
mattijs committed
213
214
215
216
217
218
219
220
221
222
223
            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;
224
            twoCells.instance() = mesh.pointsInstance();
mattijs's avatar
mattijs committed
225
226
227
228
            twoCells.write();
        }
    }

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

        if (rs.nRegions() == 1)
        {
            Info<< "    Number of regions: " << rs.nRegions() << " (OK)."
                << endl;
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
284
        }
        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"
285
286
287
288
289
290
            << setw(34) << "Surface topology";
        if (allGeometry)
        {
            Pout<< " Bounding box";
        }
        Pout<< endl;
291
292
293
294
295

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

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


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

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

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

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

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

341
                boundBox bb;   // zero-sized
342
343
344
345
346
347
348
349
350
351
352
                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>());
353
                }
354
                Pout<< ' ' << bb;
355
            }
356
            Pout<< endl;
357
358
        }

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

365
            points.instance() = mesh.pointsInstance();
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
393
            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;
}