regionSplit.C 21.9 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
OpenFOAM bot's avatar
OpenFOAM bot committed
6 7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2013 OpenFOAM Foundation
Mark Olesen's avatar
Mark Olesen committed
9
    Copyright (C) 2018-2020 OpenCFD Ltd.
10 11 12 13
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

14 15 16 17
    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.
18 19 20 21 22 23 24

    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
25
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
26 27 28 29 30 31 32 33

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

#include "regionSplit.H"
#include "cyclicPolyPatch.H"
#include "processorPolyPatch.H"
#include "globalIndex.H"
#include "syncTools.H"
34
#include "clockValue.H"
35 36 37 38 39

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

namespace Foam
{
40
    defineTypeNameAndDebug(regionSplit, 0);
41 42
}

43 44
static constexpr Foam::label UNASSIGNED = -1;
static constexpr Foam::label BLOCKED = -2;
45

mattijs's avatar
mattijs committed
46

47 48 49
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //

namespace Foam
mattijs's avatar
mattijs committed
50
{
51

52 53 54 55 56
//- The sizes of a List of containers (eg, labelHashSet)
template<class Container>
static labelList containerSizes(const UList<Container>& input)
{
    const label len = input.size();
57

58
    labelList output(len);
59

60
    for (label i = 0; i < len; ++i)
61
    {
62
        output[i] = input[i].size();
mattijs's avatar
mattijs committed
63
    }
64

65 66
    return output;
}
67

68
} // End namespace Foam
69

70

71
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
72

73 74 75 76
void Foam::regionSplit::updateFacePair
(
    const label face0,
    const label face1,
77

78 79 80 81 82 83 84 85 86 87
    labelList& faceRegion,
    DynamicList<label>& facesChanged
) const
{
    if (faceRegion[face0] == UNASSIGNED)
    {
        if (faceRegion[face1] >= 0)
        {
            faceRegion[face0] = faceRegion[face1];
            facesChanged.append(face0);
88
        }
89
        else if (faceRegion[face1] == BLOCKED)
90
        {
91 92
            // face1 blocked but not face0.
            // - illegal for coupled faces, OK for explicit connections.
93 94
        }
    }
95
    else if (faceRegion[face0] >= 0)
96
    {
97
        if (faceRegion[face1] == UNASSIGNED)
98
        {
99 100
            faceRegion[face1] = faceRegion[face0];
            facesChanged.append(face1);
101
        }
102
        else if (faceRegion[face1] == BLOCKED)
103
        {
104 105 106 107 108 109 110 111 112 113 114 115 116 117
            // face1 blocked but not face0.
            // - illegal for coupled faces, OK for explicit connections.
        }
        else if (faceRegion[face1] != faceRegion[face0])
        {
            FatalErrorInFunction
                << "Problem : coupled face " << face0
                << " on patch " << mesh().boundaryMesh().whichPatch(face0)
                << " has region " << faceRegion[face0]
                << " but coupled face " << face1
                << " has region " << faceRegion[face1] << nl
                << "Is your blocked faces specification"
                << " synchronized across coupled boundaries?" << endl
                << abort(FatalError);
118 119 120 121 122
        }
    }
}


123
void Foam::regionSplit::fillSeedMask
124
(
mattijs's avatar
mattijs committed
125
    const List<labelPair>& explicitConnections,
126 127 128 129
    const label seedCellId,
    const label markValue,
    labelList& cellRegion,
    labelList& faceRegion
130 131
) const
{
132 133
    // Seed cell
    cellRegion[seedCellId] = markValue;
134

135 136 137 138 139 140 141 142 143 144
    // Faces on seed cell
    changedFaces_.clear();
    for (const label facei : mesh().cells()[seedCellId])
    {
        if (faceRegion[facei] == UNASSIGNED)
        {
            faceRegion[facei] = markValue;
            changedFaces_.append(facei);
        }
    }
145

146 147 148 149 150
    const polyBoundaryMesh& patches = mesh().boundaryMesh();

    // Loop over changed faces. FaceCellWave in small.

    while (changedFaces_.size())
151
    {
152
        changedCells_.clear();
153

154
        for (const label facei : changedFaces_)
155
        {
156 157 158 159 160 161 162 163 164
            const label own = mesh().faceOwner()[facei];

            if (cellRegion[own] == UNASSIGNED)
            {
                cellRegion[own] = markValue;
                changedCells_.append(own);
            }

            if (mesh().isInternalFace(facei))
165
            {
166 167 168
                const label nei = mesh().faceNeighbour()[facei];

                if (cellRegion[nei] == UNASSIGNED)
169
                {
170 171
                    cellRegion[nei] = markValue;
                    changedCells_.append(nei);
172
                }
173 174 175
            }
        }

176
        if (debug & 2)
177
        {
178 179 180
            Pout<< " Changed cells / faces : "
                << changedCells_.size() << " / " << changedFaces_.size()
                << " before sync" << endl;
181 182
        }

183 184 185
        // Loop over changedCells and collect faces
        changedFaces_.clear();
        for (const label celli : changedCells_)
186
        {
187
            for (const label facei : mesh().cells()[celli])
188
            {
189 190 191 192 193
                if (faceRegion[facei] == UNASSIGNED)
                {
                    faceRegion[facei] = markValue;
                    changedFaces_.append(facei);
                }
194
            }
195 196
        }

197

198 199 200 201
        // Update locally coupled faces
        // Global connections are done later.

        for (const polyPatch& pp : patches)
202
        {
203 204 205 206 207 208 209
            if
            (
                isA<cyclicPolyPatch>(pp)
             && refCast<const cyclicPolyPatch>(pp).owner()
            )
            {
                // Transfer from neighbourPatch to here or vice versa.
210

211 212
                const cyclicPolyPatch& cycPatch =
                    refCast<const cyclicPolyPatch>(pp);
213

214
                label face0 = cycPatch.start();
215

216 217 218
                forAll(cycPatch, i)
                {
                    const label face1 = cycPatch.transformGlobalFace(face0);
219

220 221 222 223 224 225 226
                    updateFacePair
                    (
                        face0,
                        face1,
                        faceRegion,
                        changedFaces_
                    );
227

228 229 230 231
                    ++face0;
                }
            }
        }
232

233 234 235 236 237 238 239 240 241 242
        for (const labelPair& pr : explicitConnections)
        {
            updateFacePair
            (
                pr.first(),
                pr.second(),
                faceRegion,
                changedFaces_
            );
        }
243

244 245 246 247 248 249 250 251
        if (debug & 2)
        {
            Pout<< " Changed faces : "
                << changedFaces_.size()
                << " after sync" << endl;
        }
    }
}
252

253

254 255 256 257 258 259 260 261 262 263 264 265 266
Foam::label Foam::regionSplit::calcLocalRegionSplit
(
    const boolList& blockedFace,
    const List<labelPair>& explicitConnections,

    labelList& cellRegion
) const
{
    clockValue timing(debug);

    if (debug)
    {
        if (blockedFace.size())
267
        {
268 269 270
            // Check that blockedFace is synced.
            boolList syncBlockedFace(blockedFace);
            syncTools::swapFaceList(mesh(), syncBlockedFace);
271

272
            forAll(syncBlockedFace, facei)
273
            {
274
                if (syncBlockedFace[facei] != blockedFace[facei])
275
                {
276 277 278 279 280
                    FatalErrorInFunction
                        << "Face " << facei << " not synchronised. My value:"
                        << blockedFace[facei] << "  coupled value:"
                        << syncBlockedFace[facei]
                        << abort(FatalError);
281 282 283 284 285
                }
            }
        }
    }

286 287
    changedCells_.reserve(mesh_.nCells());
    changedFaces_.reserve(mesh_.nFaces());
288

289 290 291 292
    // Region per face.
    // -1 = unassigned
    // -2 = blocked
    labelList faceRegion(mesh().nFaces(), UNASSIGNED);
293

294
    if (blockedFace.size())
295
    {
296
        forAll(blockedFace, facei)
297
        {
298 299 300 301
            if (blockedFace[facei])
            {
                faceRegion[facei] = BLOCKED;
            }
302
        }
303
    }
304

305 306 307 308 309 310 311 312

    // Assign local regions
    // ~~~~~~~~~~~~~~~~~~~~

    // Start with region 0
    label nLocalRegions = 0;

    for (label seedCellId = 0; seedCellId < cellRegion.size(); ++seedCellId)
313
    {
314 315 316
        // Find next unset cell - use as seed

        for (; seedCellId < cellRegion.size(); ++seedCellId)
317
        {
318 319 320 321
            if (cellRegion[seedCellId] == UNASSIGNED)
            {
                break;
            }
322
        }
323 324

        if (seedCellId >= cellRegion.size())
325
        {
326
            break;
327 328
        }

329 330 331 332 333 334 335 336
        fillSeedMask
        (
            explicitConnections,
            seedCellId,
            nLocalRegions,
            cellRegion,
            faceRegion
        );
337

338
        ++nLocalRegions; // Next region
339 340
    }

341 342 343 344
    // Discard temporary working data
    changedCells_.clearStorage();
    changedFaces_.clearStorage();

345 346
    if (debug)
    {
347
        forAll(cellRegion, celli)
348
        {
349 350 351 352 353 354
            if (cellRegion[celli] < 0)
            {
                FatalErrorInFunction
                    << "cell:" << celli << " region:" << cellRegion[celli]
                    << abort(FatalError);
            }
355 356
        }

357
        forAll(faceRegion, facei)
358
        {
359 360 361 362 363 364
            if (faceRegion[facei] == UNASSIGNED)
            {
                FatalErrorInFunction
                    << "face:" << facei << " region:" << faceRegion[facei]
                    << abort(FatalError);
            }
365 366 367
        }
    }

Mark Olesen's avatar
Mark Olesen committed
368
    DebugInfo << "regionSplit = " << double(timing.elapsed()) << "s\n";
369

370 371
    return nLocalRegions;
}
372 373


374 375 376 377 378
Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
(
    const bool doGlobalRegions,
    const boolList& blockedFace,
    const List<labelPair>& explicitConnections,
379

380 381 382 383 384 385 386 387 388
    labelList& cellRegion
) const
{
    const label nLocalRegions = calcLocalRegionSplit
    (
        blockedFace,
        explicitConnections,
        cellRegion
    );
389

390
    if (!doGlobalRegions)
391
    {
392
        return autoPtr<globalIndex>::New(nLocalRegions);
393 394
    }

395
    return reduceRegions(nLocalRegions, blockedFace, cellRegion);
396 397 398 399 400
}


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

401 402 403 404 405
Foam::regionSplit::regionSplit
(
    const polyMesh& mesh,
    const bool doGlobalRegions
)
406
:
407
    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
mattijs's avatar
mattijs committed
408
    labelList(mesh.nCells(), -1)
409 410 411
{
    globalNumberingPtr_ = calcRegionSplit
    (
412 413 414
        doGlobalRegions,
        boolList(),         // No blockedFace
        List<labelPair>(),  // No explicitConnections
415 416 417
        *this
    );
}
418 419 420 421 422


Foam::regionSplit::regionSplit
(
    const polyMesh& mesh,
423 424
    const boolList& blockedFace,
    const bool doGlobalRegions
425 426
)
:
427
    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
mattijs's avatar
mattijs committed
428
    labelList(mesh.nCells(), -1)
429 430 431
{
    globalNumberingPtr_ = calcRegionSplit
    (
432
        doGlobalRegions,
433 434
        blockedFace,
        List<labelPair>(),  // No explicitConnections
435 436 437
        *this
    );
}
mattijs's avatar
mattijs committed
438 439 440 441 442 443


Foam::regionSplit::regionSplit
(
    const polyMesh& mesh,
    const boolList& blockedFace,
444 445
    const List<labelPair>& explicitConnections,
    const bool doGlobalRegions
mattijs's avatar
mattijs committed
446 447
)
:
448
    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
mattijs's avatar
mattijs committed
449
    labelList(mesh.nCells(), -1)
450 451 452
{
    globalNumberingPtr_ = calcRegionSplit
    (
453
        doGlobalRegions,
454 455
        blockedFace,
        explicitConnections,
456 457 458
        *this
    );
}
459 460


461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::reduceRegions
(
    const label numLocalRegions,
    const boolList& blockedFace,

    labelList& cellRegion
) const
{
    clockValue timing(debug);

    if (cellRegion.size() != mesh().nCells())
    {
        FatalErrorInFunction
            << "The cellRegion size " << cellRegion.size()
            << " is not equal to the of number of cells "
            << mesh().nCells() << endl
            << abort(FatalError);
    }


    // (numLocalRegions < 0) to signal that region information should be
    // determined ourselves. This is not really efficient, but can be useful

    const label nLocalRegions =
    (
        numLocalRegions < 0
      ? labelHashSet(cellRegion).size()
      : numLocalRegions
    );


    // Preliminary global region numbers
    const globalIndex globalRegions(nLocalRegions);


    // Lookup table of local region to global region.
    // Initially an identity mapping of the uncombined global values.

    Map<label> localToGlobal(2*nLocalRegions);
    for (const label regioni : cellRegion)
    {
        localToGlobal.insert(regioni, globalRegions.toGlobal(regioni));
    }

    // To update the localToGlobal mapping during traversal of the boundaries
    // and later when finalizing things.
    Map<label> updateLookup(2*nLocalRegions);


    // Note that we use two separate maps during the process.
    // The localToGlobal is used to map the local to global regions.
    // Merging across processors will normally make this a many->few mapping.
    // However, we may need to walk up and down processor boundaries several
    // times before all the information propagates through.
    // During these traversals, it will normally be more efficient to just
    // update the mapping without updating the cellRegion immediately.
    // Only after everything is finalized do we renumber all of the cell
    // regions.


    // Merge global regions
    // ~~~~~~~~~~~~~~~~~~~~
    // Regions across non-blocked proc patches get merged.
    // This will set merged global regions to be the min of both.
    // (this will create gaps in the global region list so they will get
    // merged later on)

    const polyBoundaryMesh& patches = mesh().boundaryMesh();

    // Buffer for swapping boundary information
533
    labelList nbrRegion(mesh().nBoundaryFaces());
534

535 536
    bool emitWarning = true;

537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601
    do
    {
        if (debug)
        {
            Pout<< nl << "-- Starting Iteration --" << endl;
        }

        updateLookup.clear();
        nbrRegion = UNASSIGNED;

        // Region information to send
        for (const polyPatch& pp : patches)
        {
            if (pp.coupled())
            {
                const labelUList& faceCells = pp.faceCells();
                SubList<label> patchNbrRegion
                (
                    nbrRegion,
                    pp.size(),
                    pp.start()-mesh().nInternalFaces()
                );

                forAll(faceCells, patchFacei)
                {
                    const label meshFacei = pp.start()+patchFacei;

                    if (!blockedFace[meshFacei])
                    {
                        // Send the most currently updated region Id
                        const label orig = cellRegion[faceCells[patchFacei]];

                        patchNbrRegion[patchFacei] = localToGlobal[orig];
                    }
                }
            }
        }
        syncTools::swapBoundaryFaceList(mesh(), nbrRegion);

        // Receive and reduce region information
        for (const polyPatch& pp : patches)
        {
            if (pp.coupled())
            {
                const labelUList& faceCells = pp.faceCells();
                SubList<label> patchNbrRegion
                (
                    nbrRegion,
                    pp.size(),
                    pp.start()-mesh().nInternalFaces()
                );

                forAll(faceCells, patchFacei)
                {
                    const label meshFacei = pp.start()+patchFacei;

                    if (!blockedFace[meshFacei])
                    {
                        // Reduction by retaining the min region id.

                        const label orig = cellRegion[faceCells[patchFacei]];

                        const label sent = localToGlobal[orig];
                        const label recv = patchNbrRegion[patchFacei];

602 603 604 605 606 607 608 609 610 611 612 613 614
                        if (recv == UNASSIGNED)
                        {
                            if (emitWarning)
                            {
                                Pout<<"Warning in regionSplit:"
                                    " received unassigned on "
                                    << pp.name() << " at patchFace "
                                    << patchFacei
                                    << ". Check synchronisation in caller"
                                    << nl;
                            }
                        }
                        else if (recv < sent)
615
                        {
616 617
                            // Record the minimum value seen

618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661
                            auto fnd = updateLookup.find(sent);
                            if (!fnd.found())
                            {
                                updateLookup.insert(sent, recv);
                            }
                            else if (recv < *fnd)
                            {
                                *fnd = recv;
                            }
                        }
                    }
                }
            }
        }


        // Note: by always using the minimum region number across the
        // processor faces, we effect a consolidation of connected regions
        // and converge to a unique number for each distinct region.


        // Update localToGlobal according to newly exchanged information

        inplaceMapValue(updateLookup, localToGlobal);

        if (debug & 2)
        {
            labelList keys(localToGlobal.sortedToc());
            labelList vals(keys.size());
            forAll(keys, i)
            {
                vals[i] = localToGlobal[keys[i]];
            }

            Pout<< "Updated local regions:" << nl
                << "old: " << flatOutput(keys) << nl
                << "new: " << flatOutput(vals) << endl;
        }
        else if (debug)
        {
            Pout<< "Updated " << localToGlobal.size()
                << " local regions" << endl;
        }

662
        emitWarning = false;
663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678
        // Continue until there are no further changes
    }
    while (returnReduce(!updateLookup.empty(), orOp<bool>()));


    //
    // We will need compact numbers locally and non-locally
    //

    // Determine the local compact numbering
    label nCompact = 0;
    {
        labelHashSet localRegion(2*localToGlobal.size());

        forAllConstIters(localToGlobal, iter)
        {
679
            const label regioni = iter.val();
680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724

            if (globalRegions.isLocal(regioni))
            {
                localRegion.insert(regioni);
            }
        }

        nCompact = localRegion.size();
    }


    // The new global numbering using compacted local regions
    auto globalCompactPtr = autoPtr<globalIndex>::New(nCompact);
    const auto& globalCompact = *globalCompactPtr;


    // Determine the following:
    // - the local compact regions (store as updateLookup)
    // - the non-local regions, ordered according to the processor on which
    //   they are local.


    // The local compaction map (updated local to compact local numbering)
    updateLookup.clear();

    labelListList sendNonLocal(Pstream::nProcs());

    {
        List<labelHashSet> nonLocal(Pstream::nProcs(), labelHashSet(0));

        // Use estimate of sizing for non-local regions
        forAll(nonLocal, proci)
        {
            if (proci != Pstream::myProcNo())
            {
                nonLocal[proci].resize
                (
                    2*((nLocalRegions-nCompact)/Pstream::nProcs())
                );
            }
        }


        forAllConstIters(localToGlobal, iter)
        {
725
            const label regioni = iter.val();
726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816

            if (globalRegions.isLocal(regioni))
            {
                updateLookup.insert
                (
                    regioni,
                    globalCompact.toGlobal(updateLookup.size())
                );
            }
            else
            {
                nonLocal[globalRegions.whichProcID(regioni)].insert(regioni);
            }
        }

        if (debug)
        {
            Pout<< " per processor nonLocal regions: "
                << flatOutput(containerSizes(nonLocal)) << endl;
        }


        // Convert to label list
        forAll(sendNonLocal, proci)
        {
            sendNonLocal[proci] = nonLocal[proci].toc();
        }
    }


    // Get the wanted region labels into recvNonLocal
    labelListList recvNonLocal;
    Pstream::exchange<labelList, label>
    (
        sendNonLocal,
        recvNonLocal
    );


    // The recvNonLocal[proci] region labels are what proci requires.
    // Transcribe into their compacted number.

    {
        labelListList sendLocal(std::move(recvNonLocal));

        for (labelList& send : sendLocal)
        {
            for (label& regioni : send)
            {
                regioni = updateLookup[regioni];
            }
        }

        // Send back (into recvNonLocal)
        Pstream::exchange<labelList, label>
        (
            sendLocal,
            recvNonLocal
        );
    }


    // Now recvNonLocal and sendNonLocal contain matched pairs with
    // sendNonLocal being the non-compact region and recvNonLocal being
    // the compact region.
    //
    // Insert these into the local compaction map.

    forAll(recvNonLocal, proci)
    {
        const labelList& send = sendNonLocal[proci];
        const labelList& recv = recvNonLocal[proci];

        forAll(send, i)
        {
            updateLookup.insert(send[i], recv[i]);
        }
    }


    // Now renumber the localToGlobal to use the final compact global values
    inplaceMapValue(updateLookup, localToGlobal);


    // Can now finally use localToGlobal to renumber cellRegion

    forAll(cellRegion, celli)
    {
        cellRegion[celli] = localToGlobal[cellRegion[celli]];
    }

Mark Olesen's avatar
Mark Olesen committed
817 818
    DebugInfo
        <<"regionSplit::reduceRegions = " << double(timing.elapsed()) << "s\n";
819 820 821 822 823

    return globalCompactPtr;
}


824
// ************************************************************************* //