regionSplit.C 12.9 KB
Newer Older
Henry's avatar
Henry committed
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
6
     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
Henry's avatar
Henry committed
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
-------------------------------------------------------------------------------
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 "regionSplit.H"
#include "cyclicPolyPatch.H"
#include "processorPolyPatch.H"
#include "globalIndex.H"
#include "syncTools.H"
31
32
#include "FaceCellWave.H"
#include "minData.H"
Henry's avatar
Henry committed
33
34
35
36
37
38
39
40
41
42
43
44

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{

defineTypeNameAndDebug(regionSplit, 0);

}

// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

45
void Foam::regionSplit::calcNonCompactRegionSplit
Henry's avatar
Henry committed
46
(
47
48
49
    const globalIndex& globalFaces,
    const boolList& blockedFace,
    const List<labelPair>& explicitConnections,
Henry's avatar
Henry committed
50

51
    labelList& cellRegion
Henry's avatar
Henry committed
52
53
) const
{
54
55
56
57
58
59
60
61
    // Field on cells and faces.
    List<minData> cellData(mesh().nCells());
    List<minData> faceData(mesh().nFaces());

    // Take over blockedFaces by seeding a negative number
    // (so is always less than the decomposition)
    label nUnblocked = 0;
    forAll(faceData, faceI)
Henry's avatar
Henry committed
62
    {
63
        if (blockedFace.size() && blockedFace[faceI])
Henry's avatar
Henry committed
64
        {
65
            faceData[faceI] = minData(-2);
Henry's avatar
Henry committed
66
        }
67
        else
Henry's avatar
Henry committed
68
        {
69
            nUnblocked++;
Henry's avatar
Henry committed
70
71
72
        }
    }

73
74
75
76
    // Seed unblocked faces
    labelList seedFaces(nUnblocked);
    List<minData> seedData(nUnblocked);
    nUnblocked = 0;
Henry's avatar
Henry committed
77
78


79
    forAll(faceData, faceI)
Henry's avatar
Henry committed
80
    {
81
        if (blockedFace.empty() || !blockedFace[faceI])
Henry's avatar
Henry committed
82
        {
83
84
85
86
            seedFaces[nUnblocked] = faceI;
            // Seed face with globally unique number
            seedData[nUnblocked] = minData(globalFaces.toGlobal(faceI));
            nUnblocked++;
Henry's avatar
Henry committed
87
88
89
90
        }
    }


91
92
93
94
95
96
97
98
99
100
101
102
    // Propagate information inwards
    FaceCellWave<minData> deltaCalc
    (
        mesh(),
        explicitConnections,
        false,  // disable walking through cyclicAMI for backwards compatibility
        seedFaces,
        seedData,
        faceData,
        cellData,
        mesh().globalData().nTotalCells()+1
    );
Henry's avatar
Henry committed
103
104


105
106
107
108
109
    // And extract
    cellRegion.setSize(mesh().nCells());
    forAll(cellRegion, cellI)
    {
        if (cellData[cellI].valid(deltaCalc.data()))
Henry's avatar
Henry committed
110
        {
111
            cellRegion[cellI] = cellData[cellI].data();
Henry's avatar
Henry committed
112
        }
113
        else
Henry's avatar
Henry committed
114
        {
115
116
            // Unvisited cell -> only possible if surrounded by blocked faces.
            // If so make up region from any of the faces
Henry's avatar
Henry committed
117
            const cell& cFaces = mesh().cells()[cellI];
118
            label faceI = cFaces[0];
Henry's avatar
Henry committed
119

120
            if (blockedFace.size() && !blockedFace[faceI])
Henry's avatar
Henry committed
121
            {
122
                FatalErrorInFunction
123
124
125
126
127
                    << "Problem: unblocked face " << faceI
                    << " at " << mesh().faceCentres()[faceI]
                    << " on unassigned cell " << cellI
                    << mesh().cellCentres()[faceI]
                    << exit(FatalError);
Henry's avatar
Henry committed
128
            }
129
            cellRegion[cellI] = globalFaces.toGlobal(faceI);
Henry's avatar
Henry committed
130
131
132
133
134
        }
    }
}


135
Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
Henry's avatar
Henry committed
136
(
137
    const bool doGlobalRegions,
Henry's avatar
Henry committed
138
139
140
141
142
143
    const boolList& blockedFace,
    const List<labelPair>& explicitConnections,

    labelList& cellRegion
) const
{
144
    // See header in regionSplit.H
Henry's avatar
Henry committed
145
146


147
    if (!doGlobalRegions)
Henry's avatar
Henry committed
148
    {
149
150
151
        // Block all parallel faces to avoid comms across
        boolList coupledOrBlockedFace(blockedFace);
        const polyBoundaryMesh& pbm = mesh().boundaryMesh();
Henry's avatar
Henry committed
152

153
        if (coupledOrBlockedFace.size())
Henry's avatar
Henry committed
154
        {
155
            forAll(pbm, patchI)
Henry's avatar
Henry committed
156
            {
157
158
159
160
161
162
163
164
165
                const polyPatch& pp = pbm[patchI];
                if (isA<processorPolyPatch>(pp))
                {
                    label faceI = pp.start();
                    forAll(pp, i)
                    {
                        coupledOrBlockedFace[faceI++] = true;
                    }
                }
Henry's avatar
Henry committed
166
167
168
            }
        }

169
170
171
        // Create dummy (local only) globalIndex
        labelList offsets(Pstream::nProcs()+1, 0);
        for (label i = Pstream::myProcNo()+1; i < offsets.size(); i++)
Henry's avatar
Henry committed
172
        {
173
            offsets[i] = mesh().nFaces();
Henry's avatar
Henry committed
174
        }
175
        const globalIndex globalRegions(offsets.xfer());
Henry's avatar
Henry committed
176

177
178
179
180
181
        // Minimise regions across connected cells
        // Note: still uses global decisions so all processors are running
        //       in lock-step, i.e. slowest determines overall time.
        //       To avoid this we could switch off Pstream::parRun.
        calcNonCompactRegionSplit
Henry's avatar
Henry committed
182
        (
183
184
            globalRegions,
            coupledOrBlockedFace,
Henry's avatar
Henry committed
185
            explicitConnections,
186
            cellRegion
Henry's avatar
Henry committed
187
188
        );

189
190
        // Compact
        Map<label> globalToCompact(mesh().nCells()/8);
Henry's avatar
Henry committed
191
192
        forAll(cellRegion, cellI)
        {
193
194
195
196
197
198
            label region = cellRegion[cellI];

            label globalRegion;

            Map<label>::const_iterator fnd = globalToCompact.find(region);
            if (fnd == globalToCompact.end())
Henry's avatar
Henry committed
199
            {
200
201
                globalRegion = globalRegions.toGlobal(globalToCompact.size());
                globalToCompact.insert(region, globalRegion);
Henry's avatar
Henry committed
202
            }
203
204
205
206
207
            else
            {
                globalRegion = fnd();
            }
            cellRegion[cellI] = globalRegion;
Henry's avatar
Henry committed
208
209
        }

210
211
212
213

        // Return globalIndex with size = localSize and all regions local
        labelList compactOffsets(Pstream::nProcs()+1, 0);
        for (label i = Pstream::myProcNo()+1; i < compactOffsets.size(); i++)
Henry's avatar
Henry committed
214
        {
215
            compactOffsets[i] = globalToCompact.size();
Henry's avatar
Henry committed
216
        }
217
218

        return autoPtr<globalIndex>(new globalIndex(compactOffsets.xfer()));
Henry's avatar
Henry committed
219
220
221
222
    }



223
224
    // Initial global region numbers
    const globalIndex globalRegions(mesh().nFaces());
Henry's avatar
Henry committed
225

226
227
    // Minimise regions across connected cells (including parallel)
    calcNonCompactRegionSplit
Henry's avatar
Henry committed
228
    (
229
        globalRegions,
Henry's avatar
Henry committed
230
231
232
233
234
235
        blockedFace,
        explicitConnections,
        cellRegion
    );


236
237
238
239
    // Now our cellRegion will have
    // - non-local regions (i.e. originating from other processors)
    // - non-compact locally originating regions
    // so we'll need to compact
Henry's avatar
Henry committed
240

241
242
    // 4a: count per originating processor the number of regions
    labelList nOriginating(Pstream::nProcs(), 0);
Henry's avatar
Henry committed
243
    {
244
        labelHashSet haveRegion(mesh().nCells()/8);
Henry's avatar
Henry committed
245

246
        forAll(cellRegion, cellI)
Henry's avatar
Henry committed
247
        {
248
            label region = cellRegion[cellI];
Henry's avatar
Henry committed
249

250
251
252
            // Count originating processor. Use isLocal as efficiency since
            // most cells are locally originating.
            if (globalRegions.isLocal(region))
Henry's avatar
Henry committed
253
            {
254
                if (haveRegion.insert(region))
Henry's avatar
Henry committed
255
                {
256
                    nOriginating[Pstream::myProcNo()]++;
Henry's avatar
Henry committed
257
258
                }
            }
259
            else
Henry's avatar
Henry committed
260
            {
261
262
                label procI = globalRegions.whichProcID(region);
                if (haveRegion.insert(region))
Henry's avatar
Henry committed
263
                {
264
                    nOriginating[procI]++;
Henry's avatar
Henry committed
265
266
267
268
269
270
271
                }
            }
        }
    }

    if (debug)
    {
272
273
        Pout<< "Counted " << nOriginating[Pstream::myProcNo()]
            << " local regions." << endl;
Henry's avatar
Henry committed
274
275
276
277
    }


    // Global numbering for compacted local regions
278
279
280
281
    autoPtr<globalIndex> globalCompactPtr
    (
        new globalIndex(nOriginating[Pstream::myProcNo()])
    );
Henry's avatar
Henry committed
282
283
284
285
286
287
288
289
290
291
    const globalIndex& globalCompact = globalCompactPtr();


    // 4b: renumber
    // Renumber into compact indices. Note that since we've already made
    // all regions global we now need a Map to store the compacting information
    // instead of a labelList - otherwise we could have used a straight
    // labelList.

    // Local compaction map
292
    Map<label> globalToCompact(2*nOriginating[Pstream::myProcNo()]);
Henry's avatar
Henry committed
293
294
295
296
    // Remote regions we want the compact number for
    List<labelHashSet> nonLocal(Pstream::nProcs());
    forAll(nonLocal, procI)
    {
297
298
299
300
        if (procI != Pstream::myProcNo())
        {
            nonLocal[procI].resize(2*nOriginating[procI]);
        }
Henry's avatar
Henry committed
301
302
303
304
305
306
307
    }

    forAll(cellRegion, cellI)
    {
        label region = cellRegion[cellI];
        if (globalRegions.isLocal(region))
        {
308
309
310
311
312
313
            // Insert new compact region (if not yet present)
            globalToCompact.insert
            (
                region,
                globalCompact.toGlobal(globalToCompact.size())
            );
Henry's avatar
Henry committed
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
        }
        else
        {
            nonLocal[globalRegions.whichProcID(region)].insert(region);
        }
    }


    // Now we have all the local regions compacted. Now we need to get the
    // non-local ones from the processors to whom they are local.
    // Convert the nonLocal (labelHashSets) to labelLists.

    labelListList sendNonLocal(Pstream::nProcs());
    forAll(sendNonLocal, procI)
    {
329
        sendNonLocal[procI] = nonLocal[procI].toc();
Henry's avatar
Henry committed
330
331
332
333
    }

    if (debug)
    {
334
        forAll(sendNonLocal, procI)
Henry's avatar
Henry committed
335
336
        {
            Pout<< "    from processor " << procI
337
                << " want " << sendNonLocal[procI].size()
Henry's avatar
Henry committed
338
339
340
341
342
343
344
345
                << " region numbers."
                << endl;
        }
        Pout<< endl;
    }


    // Get the wanted region labels into recvNonLocal
346
    labelListList recvNonLocal(Pstream::nProcs());
Henry's avatar
Henry committed
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
    labelListList sizes;
    Pstream::exchange<labelList, label>
    (
        sendNonLocal,
        recvNonLocal,
        sizes
    );

    // Now we have the wanted compact region labels that procI wants in
    // recvNonLocal[procI]. Construct corresponding list of compact
    // region labels to send back.

    labelListList sendWantedLocal(Pstream::nProcs());
    forAll(recvNonLocal, procI)
    {
        const labelList& nonLocal = recvNonLocal[procI];
        sendWantedLocal[procI].setSize(nonLocal.size());

        forAll(nonLocal, i)
        {
            sendWantedLocal[procI][i] = globalToCompact[nonLocal[i]];
        }
    }


    // Send back (into recvNonLocal)
    recvNonLocal.clear();
374
    recvNonLocal.setSize(sendWantedLocal.size());
Henry's avatar
Henry committed
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
    sizes.clear();
    Pstream::exchange<labelList, label>
    (
        sendWantedLocal,
        recvNonLocal,
        sizes
    );
    sendWantedLocal.clear();

    // Now recvNonLocal contains for every element in setNonLocal the
    // corresponding compact number. Insert these into the local compaction
    // map.

    forAll(recvNonLocal, procI)
    {
        const labelList& wantedRegions = sendNonLocal[procI];
        const labelList& compactRegions = recvNonLocal[procI];

        forAll(wantedRegions, i)
        {
            globalToCompact.insert(wantedRegions[i], compactRegions[i]);
        }
    }

    // Finally renumber the regions
    forAll(cellRegion, cellI)
    {
        cellRegion[cellI] = globalToCompact[cellRegion[cellI]];
    }

    return globalCompactPtr;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::regionSplit::regionSplit(const polyMesh& mesh, const bool doGlobalRegions)
:
    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
    labelList(mesh.nCells(), -1)
{
    globalNumberingPtr_ = calcRegionSplit
    (
        doGlobalRegions,    //do global regions
        boolList(0, false), //blockedFaces
        List<labelPair>(0), //explicitConnections,
        *this
    );
}


Foam::regionSplit::regionSplit
(
    const polyMesh& mesh,
    const boolList& blockedFace,
    const bool doGlobalRegions
)
:
    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
    labelList(mesh.nCells(), -1)
{
    globalNumberingPtr_ = calcRegionSplit
    (
        doGlobalRegions,
        blockedFace,        //blockedFaces
        List<labelPair>(0), //explicitConnections,
        *this
    );
}


Foam::regionSplit::regionSplit
(
    const polyMesh& mesh,
    const boolList& blockedFace,
    const List<labelPair>& explicitConnections,
    const bool doGlobalRegions
)
:
    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
    labelList(mesh.nCells(), -1)
{
    globalNumberingPtr_ = calcRegionSplit
    (
        doGlobalRegions,
        blockedFace,            //blockedFaces
        explicitConnections,    //explicitConnections,
        *this
    );
}


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