checkTools.C 13.4 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
6
     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
7
8
9
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
35
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "checkTools.H"
#include "polyMesh.H"
#include "globalMeshData.H"
#include "hexMatcher.H"
#include "wedgeMatcher.H"
#include "prismMatcher.H"
#include "pyrMatcher.H"
#include "tetWedgeMatcher.H"
#include "tetMatcher.H"
#include "IOmanip.H"
36
#include "pointSet.H"
37
38
39
40
41
#include "faceSet.H"
#include "cellSet.H"
#include "Time.H"
#include "surfaceWriter.H"
#include "syncTools.H"
42
43
#include "globalIndex.H"
#include "PatchTools.H"
44
#include "functionObject.H"
45

46
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
{
    Info<< "Mesh stats" << nl
        << "    points:           "
        << returnReduce(mesh.points().size(), sumOp<label>()) << nl;

    label nInternalPoints = returnReduce
    (
        mesh.nInternalPoints(),
        sumOp<label>()
    );

    if (nInternalPoints != -Pstream::nProcs())
    {
        Info<< "    internal points:  " << nInternalPoints << nl;

        if (returnReduce(mesh.nInternalPoints(), minOp<label>()) == -1)
        {
66
            WarningInFunction
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
                << "Some processors have their points sorted into internal"
                << " and external and some do not." << endl
                << "This can cause problems later on." << endl;
        }
    }

    if (allTopology && nInternalPoints != -Pstream::nProcs())
    {
        label nEdges = returnReduce(mesh.nEdges(), sumOp<label>());
        label nInternalEdges = returnReduce
        (
            mesh.nInternalEdges(),
            sumOp<label>()
        );
        label nInternal1Edges = returnReduce
        (
            mesh.nInternal1Edges(),
            sumOp<label>()
        );
        label nInternal0Edges = returnReduce
        (
            mesh.nInternal0Edges(),
            sumOp<label>()
        );

        Info<< "    edges:            " << nEdges << nl
            << "    internal edges:   " << nInternalEdges << nl
            << "    internal edges using one boundary point:   "
            << nInternal1Edges-nInternal0Edges << nl
            << "    internal edges using two boundary points:  "
            << nInternalEdges-nInternal1Edges << nl;
    }

    label nFaces = returnReduce(mesh.faces().size(), sumOp<label>());
    label nIntFaces = returnReduce(mesh.faceNeighbour().size(), sumOp<label>());
    label nCells = returnReduce(mesh.cells().size(), sumOp<label>());

    Info<< "    faces:            " << nFaces << nl
        << "    internal faces:   " << nIntFaces << nl
        << "    cells:            " << nCells << nl
        << "    faces per cell:   "
108
        << (scalar(nFaces) + scalar(nIntFaces))/max(1, nCells) << nl
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
        << "    boundary patches: " << mesh.boundaryMesh().size() << nl
        << "    point zones:      " << mesh.pointZones().size() << nl
        << "    face zones:       " << mesh.faceZones().size() << nl
        << "    cell zones:       " << mesh.cellZones().size() << nl
        << endl;

    // Construct shape recognizers
    hexMatcher hex;
    prismMatcher prism;
    wedgeMatcher wedge;
    pyrMatcher pyr;
    tetWedgeMatcher tetWedge;
    tetMatcher tet;

    // Counters for different cell types
    label nHex = 0;
    label nWedge = 0;
    label nPrism = 0;
    label nPyr = 0;
    label nTet = 0;
    label nTetWedge = 0;
    label nUnknown = 0;

    Map<label> polyhedralFaces;

134
    for (label celli = 0; celli < mesh.nCells(); celli++)
135
    {
136
        if (hex.isA(mesh, celli))
137
138
139
        {
            nHex++;
        }
140
        else if (tet.isA(mesh, celli))
141
142
143
        {
            nTet++;
        }
144
        else if (pyr.isA(mesh, celli))
145
146
147
        {
            nPyr++;
        }
148
        else if (prism.isA(mesh, celli))
149
150
151
        {
            nPrism++;
        }
152
        else if (wedge.isA(mesh, celli))
153
154
155
        {
            nWedge++;
        }
156
        else if (tetWedge.isA(mesh, celli))
157
158
159
160
161
162
        {
            nTetWedge++;
        }
        else
        {
            nUnknown++;
163
            polyhedralFaces(mesh.cells()[celli].size())++;
164
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
        }
    }

    reduce(nHex,sumOp<label>());
    reduce(nPrism,sumOp<label>());
    reduce(nWedge,sumOp<label>());
    reduce(nPyr,sumOp<label>());
    reduce(nTetWedge,sumOp<label>());
    reduce(nTet,sumOp<label>());
    reduce(nUnknown,sumOp<label>());

    Info<< "Overall number of cells of each type:" << nl
        << "    hexahedra:     " << nHex << nl
        << "    prisms:        " << nPrism << nl
        << "    wedges:        " << nWedge << nl
        << "    pyramids:      " << nPyr << nl
        << "    tet wedges:    " << nTetWedge << nl
        << "    tetrahedra:    " << nTet << nl
        << "    polyhedra:     " << nUnknown
        << endl;

    if (nUnknown > 0)
    {
        Pstream::mapCombineGather(polyhedralFaces, plusEqOp<label>());

        Info<< "    Breakdown of polyhedra by number of faces:" << nl
            << "        faces" << "   number of cells" << endl;

        const labelList sortedKeys = polyhedralFaces.sortedToc();

194
        forAll(sortedKeys, keyi)
195
        {
196
            const label nFaces = sortedKeys[keyi];
197
198
199
200
201
202
203
204
205
206
207
208

            Info<< setf(std::ios::right) << setw(13)
                << nFaces << "   " << polyhedralFaces[nFaces] << nl;
        }
    }

    Info<< endl;
}


void Foam::mergeAndWrite
(
209
    const polyMesh& mesh,
210
    const surfaceWriter& writer,
211
    const word& name,
212
    const indirectPrimitivePatch& setPatch,
213
    const fileName& outputDir
214
215
216
217
)
{
    if (Pstream::parRun())
    {
218
219
220
221
        labelList pointToGlobal;
        labelList uniqueMeshPointLabels;
        autoPtr<globalIndex> globalPoints;
        autoPtr<globalIndex> globalFaces;
222
        faceList mergedFaces;
223
224
        pointField mergedPoints;
        Foam::PatchTools::gatherAndMerge
225
        (
226
227
228
229
230
231
232
233
234
235
            mesh,
            setPatch.localFaces(),
            setPatch.meshPoints(),
            setPatch.meshPointMap(),

            pointToGlobal,
            uniqueMeshPointLabels,
            globalPoints,
            globalFaces,

236
            mergedFaces,
237
            mergedPoints
238
        );
239
240

        // Write
241
242
        if (Pstream::master())
        {
243
244
245
246
247
248
249
250
251
252
            writer.write
            (
                outputDir,
                name,
                meshedSurfRef
                (
                    mergedPoints,
                    mergedFaces
                )
            );
253
        }
254
255
256
257
258
259
    }
    else
    {
        writer.write
        (
            outputDir,
260
            name,
261
262
263
264
265
            meshedSurfRef
            (
                setPatch.localPoints(),
                setPatch.localFaces()
            )
266
267
268
269
270
        );
    }
}


271
272
273
274
275
276
277
278
279
280
281
282
283
284
void Foam::mergeAndWrite
(
    const surfaceWriter& writer,
    const faceSet& set
)
{
    const polyMesh& mesh = refCast<const polyMesh>(set.db());

    const indirectPrimitivePatch setPatch
    (
        IndirectList<face>(mesh.faces(), set.sortedToc()),
        mesh.points()
    );

285
    fileName outputDir
286
287
288
    (
        set.time().path()
      / (Pstream::parRun() ? ".." : "")
289
      / functionObject::outputPrefix
290
291
292
      / mesh.pointsInstance()
      / set.name()
    );
293
    outputDir.clean();
294
295
296
297
298

    mergeAndWrite(mesh, writer, set.name(), setPatch, outputDir);
}


299
300
301
302
303
304
305
306
307
308
309
void Foam::mergeAndWrite
(
    const surfaceWriter& writer,
    const cellSet& set
)
{
    const polyMesh& mesh = refCast<const polyMesh>(set.db());
    const polyBoundaryMesh& pbm = mesh.boundaryMesh();


    // Determine faces on outside of cellSet
310
    bitSet isInSet(mesh.nCells());
311
312
    forAllConstIter(cellSet, set, iter)
    {
313
        isInSet.set(iter.key());
314
315
316
    }


317
    boolList bndInSet(mesh.nBoundaryFaces());
318
    forAll(pbm, patchi)
319
    {
320
        const polyPatch& pp = pbm[patchi];
321
322
323
324
325
326
327
328
329
330
        const labelList& fc = pp.faceCells();
        forAll(fc, i)
        {
            bndInSet[pp.start()+i-mesh.nInternalFaces()] = isInSet[fc[i]];
        }
    }
    syncTools::swapBoundaryFaceList(mesh, bndInSet);


    DynamicList<label> outsideFaces(3*set.size());
331
    for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
332
    {
333
334
        const bool ownVal = isInSet[mesh.faceOwner()[facei]];
        const bool neiVal = isInSet[mesh.faceNeighbour()[facei]];
335
336
337

        if (ownVal != neiVal)
        {
338
            outsideFaces.append(facei);
339
340
341
342
        }
    }


343
    forAll(pbm, patchi)
344
    {
345
        const polyPatch& pp = pbm[patchi];
346
347
348
349
350
        const labelList& fc = pp.faceCells();
        if (pp.coupled())
        {
            forAll(fc, i)
            {
351
                label facei = pp.start()+i;
352

353
                const bool neiVal = bndInSet[facei-mesh.nInternalFaces()];
354
355
                if (isInSet[fc[i]] && !neiVal)
                {
356
                    outsideFaces.append(facei);
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
                }
            }
        }
        else
        {
            forAll(fc, i)
            {
                if (isInSet[fc[i]])
                {
                    outsideFaces.append(pp.start()+i);
                }
            }
        }
    }


    const indirectPrimitivePatch setPatch
    (
        IndirectList<face>(mesh.faces(), outsideFaces),
        mesh.points()
    );

379
    fileName outputDir
380
381
382
    (
        set.time().path()
      / (Pstream::parRun() ? ".." : "")
383
      / functionObject::outputPrefix
384
385
386
      / mesh.pointsInstance()
      / set.name()
    );
387
    outputDir.clean();
388

389
390
    mergeAndWrite(mesh, writer, set.name(), setPatch, outputDir);
}
391
392
393
394
395
396
397
398
399
400
401
402


void Foam::mergeAndWrite
(
    const writer<scalar>& writer,
    const pointSet& set
)
{
    const polyMesh& mesh = refCast<const polyMesh>(set.db());

    pointField mergedPts;
    labelList mergedIDs;
403
404
405

    if (Pstream::parRun())
    {
406
407
408
        // Note: we explicitly do not merge the points
        // (mesh.globalData().mergePoints etc) since this might
        // hide any synchronisation problem
409

410
411
412
413
414
415
416
417
418
        globalIndex globalNumbering(mesh.nPoints());

        mergedPts.setSize(returnReduce(set.size(), sumOp<label>()));
        mergedIDs.setSize(mergedPts.size());

        labelList setPointIDs(set.sortedToc());

        // Get renumbered local data
        pointField myPoints(mesh.points(), setPointIDs);
419
        labelList myIDs(globalNumbering.toGlobal(setPointIDs));
420

421
422
        if (Pstream::master())
        {
423
424
425
426
427
428
429
430
431
            // Insert master data first
            label pOffset = 0;
            SubList<point>(mergedPts, myPoints.size(), pOffset) = myPoints;
            SubList<label>(mergedIDs, myIDs.size(), pOffset) = myIDs;
            pOffset += myPoints.size();

            // Receive slave ones
            for (int slave=1; slave<Pstream::nProcs(); slave++)
            {
432
                IPstream fromSlave(Pstream::commsTypes::scheduled, slave);
433
434
435
436
437
438
439
440
441
442
443
444
445
446

                pointField slavePts(fromSlave);
                labelList slaveIDs(fromSlave);

                SubList<point>(mergedPts, slavePts.size(), pOffset) = slavePts;
                SubList<label>(mergedIDs, slaveIDs.size(), pOffset) = slaveIDs;
                pOffset += slaveIDs.size();
            }
        }
        else
        {
            // Construct processor stream with estimate of size. Could
            // be improved.
            OPstream toMaster
447
            (
448
                Pstream::commsTypes::scheduled,
449
450
                Pstream::masterNo(),
                myPoints.byteSize() + myIDs.byteSize()
451
            );
452
            toMaster << myPoints << myIDs;
453
        }
454
455
456
    }
    else
    {
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
        mergedIDs = set.sortedToc();
        mergedPts = pointField(mesh.points(), mergedIDs);
    }


    // Write with scalar pointID
    if (Pstream::master())
    {
        scalarField scalarPointIDs(mergedIDs.size());
        forAll(mergedIDs, i)
        {
            scalarPointIDs[i] = 1.0*mergedIDs[i];
        }

        coordSet points(set.name(), "distance", mergedPts, mag(mergedPts));

        List<const scalarField*> flds(1, &scalarPointIDs);

        wordList fldNames(1, "pointID");

        // Output e.g. pointSet p0 to
        // postProcessing/<time>/p0.vtk
479
        fileName outputDir
480
        (
481
482
            set.time().path()
          / (Pstream::parRun() ? ".." : "")
483
          / functionObject::outputPrefix
484
485
          / mesh.pointsInstance()
          // set.name()
486
        );
487
        outputDir.clean();
488
489
490
491
492
493
494
495
        mkDir(outputDir);

        fileName outputFile(outputDir/writer.getFileName(points, wordList()));
        //fileName outputFile(outputDir/set.name());

        OFstream os(outputFile);

        writer.write(points, fldNames, flds, os);
496
497
    }
}
498
499
500


// ************************************************************************* //