regionSplit.C 22.4 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 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
void Foam::regionSplit::checkBoundaryFaceSync
(
    const boolList& blockedFace
) const
{
    if (blockedFace.size())
    {
        // Check that blockedFace is synced.
        boolList syncBlockedFace(blockedFace);
        syncTools::swapFaceList(mesh(), syncBlockedFace);

        forAll(syncBlockedFace, facei)
        {
            if
            (
                blockedFace.test(facei)
             != syncBlockedFace.test(facei)
            )
            {
                FatalErrorInFunction
                    << "Face " << facei << " not synchronised. My value:"
                    << blockedFace.test(facei) << "  coupled value:"
                    << syncBlockedFace.test(facei) << nl
                    << abort(FatalError);
            }
        }
    }
}


103 104 105 106
void Foam::regionSplit::updateFacePair
(
    const label face0,
    const label face1,
107

108 109 110 111 112 113 114 115 116 117
    labelList& faceRegion,
    DynamicList<label>& facesChanged
) const
{
    if (faceRegion[face0] == UNASSIGNED)
    {
        if (faceRegion[face1] >= 0)
        {
            faceRegion[face0] = faceRegion[face1];
            facesChanged.append(face0);
118
        }
119
        else if (faceRegion[face1] == BLOCKED)
120
        {
121 122
            // face1 blocked but not face0.
            // - illegal for coupled faces, OK for explicit connections.
123 124
        }
    }
125
    else if (faceRegion[face0] >= 0)
126
    {
127
        if (faceRegion[face1] == UNASSIGNED)
128
        {
129 130
            faceRegion[face1] = faceRegion[face0];
            facesChanged.append(face1);
131
        }
132
        else if (faceRegion[face1] == BLOCKED)
133
        {
134 135 136 137 138 139 140 141 142 143 144 145 146 147
            // 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);
148 149 150 151 152
        }
    }
}


153
void Foam::regionSplit::fillSeedMask
154
(
155
    const UList<labelPair>& explicitConnections,
156 157 158 159
    const label seedCellId,
    const label markValue,
    labelList& cellRegion,
    labelList& faceRegion
160 161
) const
{
162 163
    // Seed cell
    cellRegion[seedCellId] = markValue;
164

165 166 167 168 169 170 171 172 173 174
    // Faces on seed cell
    changedFaces_.clear();
    for (const label facei : mesh().cells()[seedCellId])
    {
        if (faceRegion[facei] == UNASSIGNED)
        {
            faceRegion[facei] = markValue;
            changedFaces_.append(facei);
        }
    }
175

176 177 178 179 180
    const polyBoundaryMesh& patches = mesh().boundaryMesh();

    // Loop over changed faces. FaceCellWave in small.

    while (changedFaces_.size())
181
    {
182
        changedCells_.clear();
183

184
        for (const label facei : changedFaces_)
185
        {
186 187 188 189 190 191 192 193 194
            const label own = mesh().faceOwner()[facei];

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

            if (mesh().isInternalFace(facei))
195
            {
196 197 198
                const label nei = mesh().faceNeighbour()[facei];

                if (cellRegion[nei] == UNASSIGNED)
199
                {
200 201
                    cellRegion[nei] = markValue;
                    changedCells_.append(nei);
202
                }
203 204 205
            }
        }

206
        if (debug & 2)
207
        {
208 209 210
            Pout<< " Changed cells / faces : "
                << changedCells_.size() << " / " << changedFaces_.size()
                << " before sync" << endl;
211 212
        }

213 214 215
        // Loop over changedCells and collect faces
        changedFaces_.clear();
        for (const label celli : changedCells_)
216
        {
217
            for (const label facei : mesh().cells()[celli])
218
            {
219 220 221 222 223
                if (faceRegion[facei] == UNASSIGNED)
                {
                    faceRegion[facei] = markValue;
                    changedFaces_.append(facei);
                }
224
            }
225 226
        }

227

228 229 230 231
        // Update locally coupled faces
        // Global connections are done later.

        for (const polyPatch& pp : patches)
232
        {
233 234 235
            const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(pp);

            if (bool(cpp) && cpp->owner())
236 237
            {
                // Transfer from neighbourPatch to here or vice versa.
238
                const auto& cycPatch = *cpp;
239

240
                label face0 = cycPatch.start();
241

242 243 244
                forAll(cycPatch, i)
                {
                    const label face1 = cycPatch.transformGlobalFace(face0);
245

246 247 248 249 250 251 252
                    updateFacePair
                    (
                        face0,
                        face1,
                        faceRegion,
                        changedFaces_
                    );
253

254 255 256 257
                    ++face0;
                }
            }
        }
258

259 260 261 262 263 264 265 266 267 268
        for (const labelPair& pr : explicitConnections)
        {
            updateFacePair
            (
                pr.first(),
                pr.second(),
                faceRegion,
                changedFaces_
            );
        }
269

270 271 272 273 274 275 276 277
        if (debug & 2)
        {
            Pout<< " Changed faces : "
                << changedFaces_.size()
                << " after sync" << endl;
        }
    }
}
278

279

280
Foam::label Foam::regionSplit::localRegionSplit
281
(
282
    const UList<labelPair>& explicitConnections,
283

284 285
    labelList& cellRegion,
    labelList& faceRegion
286 287 288 289 290 291
) const
{
    clockValue timing(debug);

    changedCells_.reserve(mesh_.nCells());
    changedFaces_.reserve(mesh_.nFaces());
292

293 294 295 296 297 298 299 300

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

    // Start with region 0
    label nLocalRegions = 0;

    for (label seedCellId = 0; seedCellId < cellRegion.size(); ++seedCellId)
301
    {
302 303 304
        // Find next unset cell - use as seed

        for (; seedCellId < cellRegion.size(); ++seedCellId)
305
        {
306 307 308 309
            if (cellRegion[seedCellId] == UNASSIGNED)
            {
                break;
            }
310
        }
311 312

        if (seedCellId >= cellRegion.size())
313
        {
314
            break;
315 316
        }

317 318 319 320 321 322 323 324
        fillSeedMask
        (
            explicitConnections,
            seedCellId,
            nLocalRegions,
            cellRegion,
            faceRegion
        );
325

326
        ++nLocalRegions; // Next region
327 328
    }

329 330 331 332
    // Discard temporary working data
    changedCells_.clearStorage();
    changedFaces_.clearStorage();

333 334
    if (debug)
    {
335
        forAll(cellRegion, celli)
336
        {
337 338 339 340 341 342
            if (cellRegion[celli] < 0)
            {
                FatalErrorInFunction
                    << "cell:" << celli << " region:" << cellRegion[celli]
                    << abort(FatalError);
            }
343 344
        }

345
        forAll(faceRegion, facei)
346
        {
347 348 349 350 351 352
            if (faceRegion[facei] == UNASSIGNED)
            {
                FatalErrorInFunction
                    << "face:" << facei << " region:" << faceRegion[facei]
                    << abort(FatalError);
            }
353 354 355
        }
    }

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

358 359
    return nLocalRegions;
}
360 361


362 363
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

364 365 366 367 368
Foam::regionSplit::regionSplit
(
    const polyMesh& mesh,
    const bool doGlobalRegions
)
369
:
370
    regionSplit
371
    (
372 373
        mesh,
        bitSet(),           // No blockedFace
374
        List<labelPair>(),  // No explicitConnections
375 376 377
        doGlobalRegions
    )
{}
378 379 380 381 382


Foam::regionSplit::regionSplit
(
    const polyMesh& mesh,
383 384
    const bitSet& blockedFace,
    const List<labelPair>& explicitConnections,
385
    const bool doGlobalRegions
386 387
)
:
388
    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
389 390
    labelList(mesh.nCells(), UNASSIGNED),
    globalNumbering_()
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
    // if (debug)
    // {
    //     checkBoundaryFaceSync(blockedFace);
    // }

    labelList& cellRegion = *this;

    labelList faceRegion(mesh.nFaces(), UNASSIGNED);

    for (const label facei : blockedFace)
    {
        faceRegion[facei] = BLOCKED;
    }

    const label numLocalRegions =
        localRegionSplit(explicitConnections, cellRegion, faceRegion);

    faceRegion.clear();

    if (doGlobalRegions)
    {
        // Wrap bitset or bools
        bitSetOrBoolList hasBlockedFace(blockedFace);

        globalNumbering_ =
            reduceRegionsImpl(numLocalRegions, hasBlockedFace, cellRegion);

    }
    else
    {
        globalNumbering_ = globalIndex(numLocalRegions);
    }
424
}
mattijs's avatar
mattijs committed
425 426 427 428 429 430


Foam::regionSplit::regionSplit
(
    const polyMesh& mesh,
    const boolList& blockedFace,
431 432
    const List<labelPair>& explicitConnections,
    const bool doGlobalRegions
mattijs's avatar
mattijs committed
433 434
)
:
435
    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
436 437
    labelList(mesh.nCells(), UNASSIGNED),
    globalNumbering_()
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 468 469 470 471 472 473
    if (debug)
    {
        checkBoundaryFaceSync(blockedFace);
    }

    labelList& cellRegion = *this;

    labelList faceRegion(mesh.nFaces(), UNASSIGNED);

    forAll(blockedFace, facei)
    {
        if (blockedFace.test(facei))
        {
            faceRegion[facei] = BLOCKED;
        }
    }


    const label numLocalRegions =
        localRegionSplit(explicitConnections, cellRegion, faceRegion);

    faceRegion.clear();

    if (doGlobalRegions)
    {
        // Wrap bitset or bools
        bitSetOrBoolList hasBlockedFace(blockedFace);

        globalNumbering_ =
            reduceRegionsImpl(numLocalRegions, hasBlockedFace, cellRegion);
    }
    else
    {
        globalNumbering_ = globalIndex(numLocalRegions);
    }
474
}
475 476


477 478
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

479 480
Foam::globalIndex
Foam::regionSplit::reduceRegionsImpl
481 482
(
    const label numLocalRegions,
483
    const bitSetOrBoolList& blockedFace,
484 485 486 487 488 489 490 491 492
    labelList& cellRegion
) const
{
    clockValue timing(debug);

    if (cellRegion.size() != mesh().nCells())
    {
        FatalErrorInFunction
            << "The cellRegion size " << cellRegion.size()
493
            << " != number of cells " << mesh().nCells() << endl
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 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547
            << 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
548
    labelList nbrRegion(mesh().nBoundaryFaces());
549

550 551
    bool emitWarning = true;

552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
    do
    {
        if (debug)
        {
            Pout<< nl << "-- Starting Iteration --" << endl;
        }

        updateLookup.clear();
        nbrRegion = UNASSIGNED;

        // Region information to send
        for (const polyPatch& pp : patches)
        {
            if (pp.coupled())
            {
                SubList<label> patchNbrRegion
                (
                    nbrRegion,
                    pp.size(),
571
                    pp.offset()
572 573
                );

574
                const labelUList& faceCells = pp.faceCells();
575 576
                forAll(faceCells, patchFacei)
                {
577
                    const label celli = faceCells[patchFacei];
578 579
                    const label meshFacei = pp.start()+patchFacei;

580
                    if (!blockedFace.test(meshFacei))
581 582
                    {
                        // Send the most currently updated region Id
583
                        const label orig = cellRegion[celli];
584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600

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

        // Receive and reduce region information
        for (const polyPatch& pp : patches)
        {
            if (pp.coupled())
            {
                SubList<label> patchNbrRegion
                (
                    nbrRegion,
                    pp.size(),
601
                    pp.offset()
602 603
                );

604
                const labelUList& faceCells = pp.faceCells();
605 606
                forAll(faceCells, patchFacei)
                {
607
                    const label celli = faceCells[patchFacei];
608 609
                    const label meshFacei = pp.start()+patchFacei;

610
                    if (!blockedFace.test(meshFacei))
611 612 613
                    {
                        // Reduction by retaining the min region id.

614
                        const label orig = cellRegion[celli];
615 616 617 618

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

619 620 621 622 623 624 625 626 627 628 629 630 631
                        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)
632
                        {
633 634
                            // Record the minimum value seen

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 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678
                            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;
        }

679
        emitWarning = false;
680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695
        // 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)
        {
696
            const label regioni = iter.val();
697 698 699 700 701 702 703 704 705 706 707 708

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

        nCompact = localRegion.size();
    }


    // The new global numbering using compacted local regions
709
    globalIndex globalCompact(nCompact);
710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740


    // 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)
        {
741
            const label regioni = iter.val();
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 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832

            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
833 834
    DebugInfo
        <<"regionSplit::reduceRegions = " << double(timing.elapsed()) << "s\n";
835

836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852
    return globalCompact;
}


Foam::globalIndex
Foam::regionSplit::reduceRegions
(
    const label numLocalRegions,
    const bitSet& blockedFace,

    labelList& cellRegion
) const
{
    // Wrap bitset or bools
    bitSetOrBoolList hasBlockedFace(blockedFace);

    return reduceRegionsImpl(numLocalRegions, hasBlockedFace, cellRegion);
853 854 855
}


856
// ************************************************************************* //