globalMeshData.C 77.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 9
    Copyright (C) 2011-2017 OpenFOAM Foundation
    Copyright (C) 2015-2019 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 "globalMeshData.H"
#include "Pstream.H"
#include "PstreamCombineReduceOps.H"
#include "processorPolyPatch.H"
#include "globalPoints.H"
34 35
#include "polyMesh.H"
#include "mapDistribute.H"
36 37
#include "labelIOList.H"
#include "mergePoints.H"
38
#include "globalIndexAndTransform.H"
39 40 41

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

42 43 44
namespace Foam
{
defineTypeNameAndDebug(globalMeshData, 0);
45

46
const scalar globalMeshData::matchTol_ = 1e-8;
47

48 49 50 51 52 53 54 55 56 57 58 59
template<>
class minEqOp<labelPair>
{
public:
    void operator()(labelPair& x, const labelPair& y) const
    {
        x[0] = min(x[0], y[0]);
        x[1] = min(x[1], y[1]);
    }
};
}

60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76

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

void Foam::globalMeshData::initProcAddr()
{
    processorPatchIndices_.setSize(mesh_.boundaryMesh().size());
    processorPatchIndices_ = -1;

    processorPatchNeighbours_.setSize(mesh_.boundaryMesh().size());
    processorPatchNeighbours_ = -1;

    // Construct processor patch indexing. processorPatchNeighbours_ only
    // set if running in parallel!
    processorPatches_.setSize(mesh_.boundaryMesh().size());

    label nNeighbours = 0;

77
    forAll(mesh_.boundaryMesh(), patchi)
78 79 80 81 82 83 84 85 86 87 88 89
    {
        if (isA<processorPolyPatch>(mesh_.boundaryMesh()[patchi]))
        {
            processorPatches_[nNeighbours] = patchi;
            processorPatchIndices_[patchi] = nNeighbours++;
        }
    }
    processorPatches_.setSize(nNeighbours);


    if (Pstream::parRun())
    {
90
        PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
91

92
        // Send indices of my processor patches to my neighbours
93
        for (const label patchi : processorPatches_)
94
        {
95
            UOPstream toNeighbour
96 97 98 99
            (
                refCast<const processorPolyPatch>
                (
                    mesh_.boundaryMesh()[patchi]
100 101
                ).neighbProcNo(),
                pBufs
102 103 104 105 106
            );

            toNeighbour << processorPatchIndices_[patchi];
        }

107 108
        pBufs.finishedSends();

109
        for (const label patchi : processorPatches_)
110
        {
111
            UIPstream fromNeighbour
112 113 114 115
            (
                refCast<const processorPolyPatch>
                (
                    mesh_.boundaryMesh()[patchi]
116 117
                ).neighbProcNo(),
                pBufs
118
            );
119

120 121 122 123 124 125
            fromNeighbour >> processorPatchNeighbours_[patchi];
        }
    }
}


126 127 128 129 130 131 132 133 134
void Foam::globalMeshData::calcSharedPoints() const
{
    if
    (
        nGlobalPoints_ != -1
     || sharedPointLabelsPtr_.valid()
     || sharedPointAddrPtr_.valid()
    )
    {
135
        FatalErrorInFunction
136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 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 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
            << "Shared point addressing already done" << abort(FatalError);
    }

    // Calculate all shared points (exclude points that are only
    // on two coupled patches). This does all the hard work.
    globalPoints parallelPoints(mesh_, false, true);

    // Count the number of master points
    label nMaster = 0;
    forAll(parallelPoints.pointPoints(), i)
    {
        const labelList& pPoints = parallelPoints.pointPoints()[i];
        const labelList& transPPoints =
            parallelPoints.transformedPointPoints()[i];

        if (pPoints.size()+transPPoints.size() > 0)
        {
            nMaster++;
        }
    }

    // Allocate global numbers
    globalIndex masterNumbering(nMaster);

    nGlobalPoints_ = masterNumbering.size();


    // Push master number to slaves
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // 1. Fill master and slave slots
    nMaster = 0;
    labelList master(parallelPoints.map().constructSize(), -1);
    forAll(parallelPoints.pointPoints(), i)
    {
        const labelList& pPoints = parallelPoints.pointPoints()[i];
        const labelList& transPPoints =
            parallelPoints.transformedPointPoints()[i];

        if (pPoints.size()+transPPoints.size() > 0)
        {
            master[i] = masterNumbering.toGlobal(nMaster);
            forAll(pPoints, j)
            {
                master[pPoints[j]] = master[i];
            }
            forAll(transPPoints, j)
            {
                master[transPPoints[j]] = master[i];
            }
            nMaster++;
        }
    }


    // 2. Push slave slots back to local storage on originating processor
    // For all the four types of points:
    // - local master : already set
    // - local transformed slave point : the reverse transform at
    //   reverseDistribute will have copied it back to its originating local
    //   point
    // - remote untransformed slave point : sent back to originating processor
    // - remote transformed slave point : the reverse transform will
    //   copy it back into the remote slot which then gets sent back to
    //   originating processor

    parallelPoints.map().reverseDistribute
    (
        parallelPoints.map().constructSize(),
        master
    );


    // Collect all points that are a master or refer to a master.
    nMaster = 0;
    forAll(parallelPoints.pointPoints(), i)
    {
        if (master[i] != -1)
        {
            nMaster++;
        }
    }

    sharedPointLabelsPtr_.reset(new labelList(nMaster));
    labelList& sharedPointLabels = sharedPointLabelsPtr_();
    sharedPointAddrPtr_.reset(new labelList(nMaster));
    labelList& sharedPointAddr = sharedPointAddrPtr_();
    nMaster = 0;

    forAll(parallelPoints.pointPoints(), i)
    {
        if (master[i] != -1)
        {
            // I am master or slave
            sharedPointLabels[nMaster] = i;
            sharedPointAddr[nMaster] = master[i];
            nMaster++;
        }
    }

    if (debug)
    {
        Pout<< "globalMeshData : nGlobalPoints_:" << nGlobalPoints_ << nl
            << "globalMeshData : sharedPointLabels_:"
            << sharedPointLabelsPtr_().size() << nl
            << "globalMeshData : sharedPointAddr_:"
            << sharedPointAddrPtr_().size() << endl;
    }
}


246 247
void Foam::globalMeshData::countSharedEdges
(
248 249
    const EdgeMap<labelList>& procSharedEdges,
    EdgeMap<label>& globalShared,
250 251 252 253
    label& sharedEdgeI
)
{
    // Count occurrences of procSharedEdges in global shared edges table.
254
    forAllConstIters(procSharedEdges, iter)
255 256 257
    {
        const edge& e = iter.key();

258
        auto globalFnd = globalShared.find(e);
259

260 261 262 263
        if (globalFnd.found())
        {
            if (globalFnd() == -1)
            {
264
                // Second time occurrence of this edge.
265 266 267 268 269
                // Assign proper edge label.
                globalFnd() = sharedEdgeI++;
            }
        }
        else
270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289
        {
            // First time occurrence of this edge. Check how many we are adding.
            if (iter().size() == 1)
            {
                // Only one edge. Mark with special value.
                globalShared.insert(e, -1);
            }
            else
            {
                // Edge used more than once (even by local shared edges alone)
                // so allocate proper shared edge label.
                globalShared.insert(e, sharedEdgeI++);
            }
        }
    }
}


void Foam::globalMeshData::calcSharedEdges() const
{
290 291 292 293 294
    // Shared edges are shared between multiple processors. By their nature both
    // of their endpoints are shared points. (but not all edges using two shared
    // points are shared edges! There might e.g. be an edge between two
    // unrelated clusters of shared points)

295 296 297 298 299 300
    if
    (
        nGlobalEdges_ != -1
     || sharedEdgeLabelsPtr_.valid()
     || sharedEdgeAddrPtr_.valid()
    )
301
    {
302
        FatalErrorInFunction
303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320
            << "Shared edge addressing already done" << abort(FatalError);
    }


    const labelList& sharedPtAddr = sharedPointAddr();
    const labelList& sharedPtLabels = sharedPointLabels();

    // Since don't want to construct pointEdges for whole mesh create
    // Map for all shared points.
    Map<label> meshToShared(2*sharedPtLabels.size());
    forAll(sharedPtLabels, i)
    {
        meshToShared.insert(sharedPtLabels[i], i);
    }

    // Find edges using shared points. Store correspondence to local edge
    // numbering. Note that multiple local edges can have the same shared
    // points! (for cyclics or separated processor patches)
321
    EdgeMap<labelList> localShared(2*sharedPtAddr.size());
322 323 324 325 326 327 328

    const edgeList& edges = mesh_.edges();

    forAll(edges, edgeI)
    {
        const edge& e = edges[edgeI];

329
        const auto e0Fnd = meshToShared.cfind(e[0]);
330

331
        if (e0Fnd.found())
332
        {
333
            const auto e1Fnd = meshToShared.cfind(e[1]);
334

335
            if (e1Fnd.found())
336 337 338 339 340 341 342 343 344 345 346
            {
                // Found edge which uses shared points. Probably shared.

                // Construct the edge in shared points (or rather global indices
                // of the shared points)
                edge sharedEdge
                (
                    sharedPtAddr[e0Fnd()],
                    sharedPtAddr[e1Fnd()]
                );

347
                auto iter = localShared.find(sharedEdge);
348

349
                if (!iter.found())
350 351 352 353 354 355 356 357 358
                {
                    // First occurrence of this point combination. Store.
                    localShared.insert(sharedEdge, labelList(1, edgeI));
                }
                else
                {
                    // Add this edge to list of edge labels.
                    labelList& edgeLabels = iter();

359
                    const label sz = edgeLabels.size();
360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376
                    edgeLabels.setSize(sz+1);
                    edgeLabels[sz] = edgeI;
                }
            }
        }
    }


    // Now we have a table on every processors which gives its edges which use
    // shared points. Send this all to the master and have it allocate
    // global edge numbers for it. But only allocate a global edge number for
    // edge if it is used more than once!
    // Note that we are now sending the whole localShared to the master whereas
    // we only need the local count (i.e. the number of times a global edge is
    // used). But then this only gets done once so not too bothered about the
    // extra global communication.

377
    EdgeMap<label> globalShared(nGlobalPoints());
378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401

    if (Pstream::master())
    {
        label sharedEdgeI = 0;

        // Merge my shared edges into the global list
        if (debug)
        {
            Pout<< "globalMeshData::calcSharedEdges : Merging in from proc0 : "
                << localShared.size() << endl;
        }
        countSharedEdges(localShared, globalShared, sharedEdgeI);

        // Receive data from slaves and insert
        if (Pstream::parRun())
        {
            for
            (
                int slave=Pstream::firstSlave();
                slave<=Pstream::lastSlave();
                slave++
            )
            {
                // Receive the edges using shared points from the slave.
402
                IPstream fromSlave(Pstream::commsTypes::blocking, slave);
403
                EdgeMap<labelList> procSharedEdges(fromSlave);
404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419

                if (debug)
                {
                    Pout<< "globalMeshData::calcSharedEdges : "
                        << "Merging in from proc"
                        << Foam::name(slave) << " : " << procSharedEdges.size()
                        << endl;
                }
                countSharedEdges(procSharedEdges, globalShared, sharedEdgeI);
            }
        }

        // Now our globalShared should have some edges with -1 as edge label
        // These were only used once so are not proper shared edges.
        // Remove them.
        {
420
            EdgeMap<label> oldSharedEdges(globalShared);
421 422 423

            globalShared.clear();

424
            forAllConstIters(oldSharedEdges, iter)
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
            {
                if (iter() != -1)
                {
                    globalShared.insert(iter.key(), iter());
                }
            }
            if (debug)
            {
                Pout<< "globalMeshData::calcSharedEdges : Filtered "
                    << oldSharedEdges.size()
                    << " down to " << globalShared.size() << endl;
            }
        }


        // Send back to slaves.
        if (Pstream::parRun())
        {
            for
            (
                int slave=Pstream::firstSlave();
                slave<=Pstream::lastSlave();
                slave++
            )
            {
                // Receive the edges using shared points from the slave.
451
                OPstream toSlave(Pstream::commsTypes::blocking, slave);
452 453 454 455 456 457 458 459
                toSlave << globalShared;
            }
        }
    }
    else
    {
        // Send local edges to master
        {
460 461 462 463 464
            OPstream toMaster
            (
                Pstream::commsTypes::blocking,
                Pstream::masterNo()
            );
465 466 467 468
            toMaster << localShared;
        }
        // Receive merged edges from master.
        {
469 470 471 472 473
            IPstream fromMaster
            (
                Pstream::commsTypes::blocking,
                Pstream::masterNo()
            );
474 475 476 477 478 479 480 481 482 483 484 485
            fromMaster >> globalShared;
        }
    }

    // Now use the global shared edges list (globalShared) to classify my local
    // ones (localShared)

    nGlobalEdges_ = globalShared.size();

    DynamicList<label> dynSharedEdgeLabels(globalShared.size());
    DynamicList<label> dynSharedEdgeAddr(globalShared.size());

486
    forAllConstIters(localShared, iter)
487 488 489
    {
        const edge& e = iter.key();

490
        const auto edgeFnd = globalShared.cfind(e);
491

492
        if (edgeFnd.found())
493 494 495 496 497
        {
            // My local edge is indeed a shared one. Go through all local edge
            // labels with this point combination.
            const labelList& edgeLabels = iter();

498
            for (const label edgei : edgeLabels)
499 500
            {
                // Store label of local mesh edge
501
                dynSharedEdgeLabels.append(edgei);
502 503 504 505 506 507

                // Store label of shared edge
                dynSharedEdgeAddr.append(edgeFnd());
            }
        }
    }
508

509

510 511
    sharedEdgeLabelsPtr_.reset(new labelList());
    labelList& sharedEdgeLabels = sharedEdgeLabelsPtr_();
512 513
    sharedEdgeLabels.transfer(dynSharedEdgeLabels);

514 515
    sharedEdgeAddrPtr_.reset(new labelList());
    labelList& sharedEdgeAddr = sharedEdgeAddrPtr_();
516 517 518 519 520 521 522 523 524 525 526 527 528
    sharedEdgeAddr.transfer(dynSharedEdgeAddr);

    if (debug)
    {
        Pout<< "globalMeshData : nGlobalEdges_:" << nGlobalEdges_ << nl
            << "globalMeshData : sharedEdgeLabels:" << sharedEdgeLabels.size()
            << nl
            << "globalMeshData : sharedEdgeAddr:" << sharedEdgeAddr.size()
            << endl;
    }
}


529
void Foam::globalMeshData::calcGlobalPointSlaves() const
530
{
531
    if (debug)
532
    {
533 534 535
        Pout<< "globalMeshData::calcGlobalPointSlaves() :"
            << " calculating coupled master to slave point addressing."
            << endl;
536 537
    }

538 539
    // Calculate connected points for master points.
    globalPoints globalData(mesh_, coupledPatch(), true, true);
540

541
    globalPointSlavesPtr_.reset
542
    (
543 544
        new labelListList
        (
545
            std::move(globalData.pointPoints())
546
        )
547
    );
548
    globalPointTransformedSlavesPtr_.reset
549
    (
550
        new labelListList
551
        (
552
            std::move(globalData.transformedPointPoints())
553
        )
554
    );
555

556
    globalPointSlavesMapPtr_.reset
557 558 559
    (
        new mapDistribute
        (
560
            std::move(globalData.map())
561 562 563 564 565
        )
    );
}


566 567 568 569
void Foam::globalMeshData::calcPointConnectivity
(
    List<labelPairList>& allPointConnectivity
) const
570
{
571 572 573 574 575 576 577
    const globalIndexAndTransform& transforms = globalTransforms();
    const labelListList& slaves = globalPointSlaves();
    const labelListList& transformedSlaves = globalPointTransformedSlaves();


    // Create field with my local data
    labelPairList myData(globalPointSlavesMap().constructSize());
578
    forAll(slaves, pointi)
579
    {
580
        myData[pointi] = transforms.encode
581 582
        (
            Pstream::myProcNo(),
583
            pointi,
584 585
            transforms.nullTransformIndex()
        );
586
    }
587
    // Send to master
588
    globalPointSlavesMap().distribute(myData);
589 590


591 592
    // String of connected points with their transform
    allPointConnectivity.setSize(globalPointSlavesMap().constructSize());
593 594 595 596
    allPointConnectivity = labelPairList(0);

    // Pass1: do the master points since these also update local slaves
    //        (e.g. from local cyclics)
597
    forAll(slaves, pointi)
598 599
    {
        // Reconstruct string of connected points
600 601
        const labelList& pSlaves = slaves[pointi];
        const labelList& pTransformSlaves = transformedSlaves[pointi];
602

603
        if (pSlaves.size()+pTransformSlaves.size())
604
        {
605
            labelPairList& pConnectivity = allPointConnectivity[pointi];
606

607 608 609 610
            pConnectivity.setSize(1+pSlaves.size()+pTransformSlaves.size());
            label connI = 0;

            // Add myself
611
            pConnectivity[connI++] = myData[pointi];
612 613 614 615 616 617 618 619 620 621 622 623 624 625 626
            // Add untransformed points
            forAll(pSlaves, i)
            {
                pConnectivity[connI++] = myData[pSlaves[i]];
            }
            // Add transformed points.
            forAll(pTransformSlaves, i)
            {
                // Get transform from index
                label transformI = globalPointSlavesMap().whichTransform
                (
                    pTransformSlaves[i]
                );
                // Add transform to connectivity
                const labelPair& n = myData[pTransformSlaves[i]];
627 628 629
                label proci = transforms.processor(n);
                label index = transforms.index(n);
                pConnectivity[connI++] = transforms.encode
630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645
                (
                    proci,
                    index,
                    transformI
                );
            }

            // Put back in slots
            forAll(pSlaves, i)
            {
                allPointConnectivity[pSlaves[i]] = pConnectivity;
            }
            forAll(pTransformSlaves, i)
            {
                allPointConnectivity[pTransformSlaves[i]] = pConnectivity;
            }
646
        }
647 648 649 650
    }


    // Pass2: see if anything is still unset (should not be the case)
651
    forAll(slaves, pointi)
652
    {
653
        labelPairList& pConnectivity = allPointConnectivity[pointi];
654 655

        if (pConnectivity.size() == 0)
656
        {
657
            pConnectivity.setSize(1, myData[pointi]);
658 659
        }
    }
660 661


662 663
    globalPointSlavesMap().reverseDistribute
    (
664
        slaves.size(),
665 666 667
        allPointConnectivity
    );
}
668 669


670 671 672 673 674 675 676
void Foam::globalMeshData::calcGlobalPointEdges
(
    labelListList& globalPointEdges,
    List<labelPairList>& globalPointPoints
) const
{
    const edgeList& edges = coupledPatch().edges();
677
    const labelListList& pointEdges = coupledPatch().pointEdges();
678
    const globalIndex& globalEdgeNumbers = globalEdgeNumbering();
679 680
    const labelListList& slaves = globalPointSlaves();
    const labelListList& transformedSlaves = globalPointTransformedSlaves();
681 682
    const globalIndexAndTransform& transforms = globalTransforms();

683 684 685 686

    // Create local version
    globalPointEdges.setSize(globalPointSlavesMap().constructSize());
    globalPointPoints.setSize(globalPointSlavesMap().constructSize());
687
    forAll(pointEdges, pointi)
688
    {
689
        const labelList& pEdges = pointEdges[pointi];
690
        globalPointEdges[pointi] = globalEdgeNumbers.toGlobal(pEdges);
691

692
        labelPairList& globalPPoints = globalPointPoints[pointi];
693 694 695
        globalPPoints.setSize(pEdges.size());
        forAll(pEdges, i)
        {
696
            label otherPointi = edges[pEdges[i]].otherVertex(pointi);
697
            globalPPoints[i] = transforms.encode
698 699
            (
                Pstream::myProcNo(),
700
                otherPointi,
701
                transforms.nullTransformIndex()
702 703
            );
        }
704 705
    }

706 707
    // Pull slave data to master. Dummy transform.
    globalPointSlavesMap().distribute(globalPointEdges);
708 709
    globalPointSlavesMap().distribute(globalPointPoints);
    // Add all pointEdges
710
    forAll(slaves, pointi)
711
    {
712 713
        const labelList& pSlaves = slaves[pointi];
        const labelList& pTransformSlaves = transformedSlaves[pointi];
714

715 716 717 718 719 720 721 722 723
        label n = 0;
        forAll(pSlaves, i)
        {
            n += globalPointEdges[pSlaves[i]].size();
        }
        forAll(pTransformSlaves, i)
        {
            n += globalPointEdges[pTransformSlaves[i]].size();
        }
724

725 726
        // Add all the point edges of the slaves to those of the (master) point
        {
727
            labelList& globalPEdges = globalPointEdges[pointi];
728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746
            label sz = globalPEdges.size();
            globalPEdges.setSize(sz+n);
            forAll(pSlaves, i)
            {
                const labelList& otherData = globalPointEdges[pSlaves[i]];
                forAll(otherData, j)
                {
                    globalPEdges[sz++] = otherData[j];
                }
            }
            forAll(pTransformSlaves, i)
            {
                const labelList& otherData =
                    globalPointEdges[pTransformSlaves[i]];
                forAll(otherData, j)
                {
                    globalPEdges[sz++] = otherData[j];
                }
            }
747

748 749 750 751 752 753 754 755 756 757
            // Put back in slots
            forAll(pSlaves, i)
            {
                globalPointEdges[pSlaves[i]] = globalPEdges;
            }
            forAll(pTransformSlaves, i)
            {
                globalPointEdges[pTransformSlaves[i]] = globalPEdges;
            }
        }
758

759

760
        // Same for corresponding pointPoints
761
        {
762
            labelPairList& globalPPoints = globalPointPoints[pointi];
763 764
            label sz = globalPPoints.size();
            globalPPoints.setSize(sz + n);
765

766 767
            // Add untransformed points
            forAll(pSlaves, i)
768
            {
769 770
                const labelPairList& otherData = globalPointPoints[pSlaves[i]];
                forAll(otherData, j)
771
                {
772
                    globalPPoints[sz++] = otherData[j];
773 774
                }
            }
775 776
            // Add transformed points.
            forAll(pTransformSlaves, i)
777
            {
778 779 780 781 782
                // Get transform from index
                label transformI = globalPointSlavesMap().whichTransform
                (
                    pTransformSlaves[i]
                );
783

784 785 786 787 788 789
                const labelPairList& otherData =
                    globalPointPoints[pTransformSlaves[i]];
                forAll(otherData, j)
                {
                    // Add transform to connectivity
                    const labelPair& n = otherData[j];
790 791 792
                    label proci = transforms.processor(n);
                    label index = transforms.index(n);
                    globalPPoints[sz++] = transforms.encode
793
                    (
794
                        proci,
795 796 797
                        index,
                        transformI
                    );
798 799
                }
            }
800 801 802 803 804 805 806 807 808 809

            // Put back in slots
            forAll(pSlaves, i)
            {
                globalPointPoints[pSlaves[i]] = globalPPoints;
            }
            forAll(pTransformSlaves, i)
            {
                globalPointPoints[pTransformSlaves[i]] = globalPPoints;
            }
810
        }
811 812 813 814
    }
    // Push back
    globalPointSlavesMap().reverseDistribute
    (
815
        slaves.size(),
816 817 818 819 820
        globalPointEdges
    );
    // Push back
    globalPointSlavesMap().reverseDistribute
    (
821
        slaves.size(),
822 823 824
        globalPointPoints
    );
}
825 826


827 828 829 830 831 832 833
Foam::label Foam::globalMeshData::findTransform
(
    const labelPairList& info,
    const labelPair& remotePoint,
    const label localPoint
) const
{
834 835 836 837
    const globalIndexAndTransform& transforms = globalTransforms();

    const label remoteProci = transforms.processor(remotePoint);
    const label remoteIndex = transforms.index(remotePoint);
838 839 840 841 842

    label remoteTransformI = -1;
    label localTransformI = -1;
    forAll(info, i)
    {
843 844 845
        label proci = transforms.processor(info[i]);
        label pointi = transforms.index(info[i]);
        label transformI = transforms.transformIndex(info[i]);
846

847
        if (proci == Pstream::myProcNo() && pointi == localPoint)
848 849 850 851 852 853
        {
            localTransformI = transformI;
            //Pout<< "For local :" << localPoint
            //    << " found transform:" << localTransformI
            //    << endl;
        }
854
        if (proci == remoteProci && pointi == remoteIndex)
855
        {
856 857 858 859 860 861 862
            remoteTransformI = transformI;
            //Pout<< "For remote:" << remotePoint
            //    << " found transform:" << remoteTransformI
            //    << " at index:" << i
            //    << endl;
        }
    }
863

864 865
    if (remoteTransformI == -1 || localTransformI == -1)
    {
866
        FatalErrorInFunction
867
            << "Problem. Cannot find " << remotePoint
868 869 870
            << " or " << localPoint  << " "
            << coupledPatch().localPoints()[localPoint]
            << " in " << info
871 872 873 874 875
            << endl
            << "remoteTransformI:" << remoteTransformI << endl
            << "localTransformI:" << localTransformI
            << abort(FatalError);
    }
876

877
    return transforms.subtractTransformIndex
878 879 880 881 882
    (
        remoteTransformI,
        localTransformI
    );
}
883

884 885 886 887 888 889 890 891 892 893 894

void Foam::globalMeshData::calcGlobalEdgeSlaves() const
{
    if (debug)
    {
        Pout<< "globalMeshData::calcGlobalEdgeSlaves() :"
            << " calculating coupled master to slave edge addressing." << endl;
    }

    const edgeList& edges = coupledPatch().edges();
    const globalIndex& globalEdgeNumbers = globalEdgeNumbering();
895
    const globalIndexAndTransform& transforms = globalTransforms();
896 897 898


    // The whole problem with deducting edge-connectivity from
Andrew Heather's avatar
Andrew Heather committed
899
    // point-connectivity is that one of the endpoints might be
900 901 902 903 904
    // a local master but the other endpoint might not. So we first
    // need to make sure that all points know about connectivity and
    // the transformations.


laurence's avatar
laurence committed
905
    // 1. collect point connectivity - basically recreating globalPoints output.
906
    // All points will now have a string of coupled points. The transforms are
907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943
    // in respect to the master.
    List<labelPairList> allPointConnectivity;
    calcPointConnectivity(allPointConnectivity);


    // 2. Get all pointEdges and pointPoints
    // Coupled point to global coupled edges and corresponding endpoint.
    labelListList globalPointEdges;
    List<labelPairList> globalPointPoints;
    calcGlobalPointEdges(globalPointEdges, globalPointPoints);


    // 3. Now all points have
    //      - all the connected points with original transform
    //      - all the connected global edges

    // Now all we need to do is go through all the edges and check
    // both endpoints. If there is a edge between the two which is
    // produced by transforming both points in the same way it is a shared
    // edge.

    // Collect strings of connected edges.
    List<labelPairList> allEdgeConnectivity(edges.size());

    forAll(edges, edgeI)
    {
        const edge& e = edges[edgeI];
        const labelList& pEdges0 = globalPointEdges[e[0]];
        const labelPairList& pPoints0 = globalPointPoints[e[0]];
        const labelList& pEdges1 = globalPointEdges[e[1]];
        const labelPairList& pPoints1 = globalPointPoints[e[1]];

        // Most edges will be size 2
        DynamicList<labelPair> eEdges(2);
        // Append myself.
        eEdges.append
        (
944
            transforms.encode
945 946 947
            (
                Pstream::myProcNo(),
                edgeI,
948
                transforms.nullTransformIndex()
949 950 951 952 953 954
            )
        );

        forAll(pEdges0, i)
        {
            forAll(pEdges1, j)
955
            {
956 957 958 959 960
                if
                (
                    pEdges0[i] == pEdges1[j]
                 && pEdges0[i] != globalEdgeNumbers.toGlobal(edgeI)
                )
961
                {
962 963 964 965 966
                    // Found a shared edge. Now check if the endpoints
                    // go through the same transformation.
                    // Local: e[0]    remote:pPoints1[j]
                    // Local: e[1]    remote:pPoints0[i]

967

968 969 970 971
                    // Find difference in transforms to go from point on remote
                    // edge (pPoints1[j]) to this point.

                    label transform0 = findTransform
972
                    (
973 974 975 976 977 978 979 980 981
                        allPointConnectivity[e[0]],
                        pPoints1[j],
                        e[0]
                    );
                    label transform1 = findTransform
                    (
                        allPointConnectivity[e[1]],
                        pPoints0[i],
                        e[1]
982 983
                    );

984 985
                    if (transform0 == transform1)
                    {
986
                        label proci = globalEdgeNumbers.whichProcID(pEdges0[i]);
987
                        eEdges.append
988
                        (
989
                            transforms.encode
990
                            (
991 992
                                proci,
                                globalEdgeNumbers.toLocal(proci, pEdges0[i]),
993
                                transform0
994 995
                            )
                        );
996 997 998 999 1000 1001
                    }
                }
            }
        }

        allEdgeConnectivity[edgeI].transfer(eEdges);
1002 1003 1004 1005 1006
        sort
        (
            allEdgeConnectivity[edgeI],
            globalIndexAndTransform::less(transforms)
        );
1007 1008 1009 1010 1011
    }

    // We now have - in allEdgeConnectivity - a list of edges which are shared
    // between multiple processors. Filter into non-transformed and transformed
    // connections.
1012

1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026
    globalEdgeSlavesPtr_.reset(new labelListList(edges.size()));
    labelListList& globalEdgeSlaves = globalEdgeSlavesPtr_();
    List<labelPairList> transformedEdges(edges.size());
    forAll(allEdgeConnectivity, edgeI)
    {
        const labelPairList& edgeInfo = allEdgeConnectivity[edgeI];
        if (edgeInfo.size() >= 2)
        {
            const labelPair& masterInfo = edgeInfo[0];

            // Check if master edge (= first element (since sorted)) is me.
            if
            (
                (
1027
                    transforms.processor(masterInfo)
1028 1029
                 == Pstream::myProcNo()
                )
1030
             && (transforms.index(masterInfo) == edgeI)
1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045
            )
            {
                // Sort into transformed and untransformed
                labelList& eEdges = globalEdgeSlaves[edgeI];
                eEdges.setSize(edgeInfo.size()-1);

                labelPairList& trafoEEdges = transformedEdges[edgeI];
                trafoEEdges.setSize(edgeInfo.size()-1);

                label nonTransformI = 0;
                label transformI = 0;

                for (label i = 1; i < edgeInfo.size(); i++)
                {
                    const labelPair& info = edgeInfo[i];
1046 1047 1048
                    label proci = transforms.processor(info);
                    label index = transforms.index(info);
                    label transform = transforms.transformIndex
1049 1050 1051 1052
                    (
                        info
                    );

1053
                    if (transform == transforms.nullTransformIndex())
1054 1055 1056
                    {
                        eEdges[nonTransformI++] = globalEdgeNumbers.toGlobal
                        (
1057
                            proci,
1058 1059 1060 1061 1062 1063
                            index
                        );
                    }
                    else
                    {
                        trafoEEdges[transformI++] = info;
1064 1065
                    }
                }
1066 1067 1068

                eEdges.setSize(nonTransformI);
                trafoEEdges.setSize(transformI);
1069 1070 1071 1072 1073
            }
        }
    }


1074 1075
    // Construct map
    globalEdgeTransformedSlavesPtr_.reset(new labelListList());
1076

1077
    List<Map<label>> compactMap(Pstream::nProcs());
1078
    globalEdgeSlavesMapPtr_.reset
1079 1080 1081
    (
        new mapDistribute
        (
1082 1083
            globalEdgeNumbers,
            globalEdgeSlaves,
mattijs's avatar
mattijs committed
1084

1085
            transforms,
1086 1087
            transformedEdges,
            globalEdgeTransformedSlavesPtr_(),
mattijs's avatar
mattijs committed
1088

1089 1090
            compactMap
        )
mattijs's avatar
mattijs committed
1091 1092 1093 1094 1095
    );


    if (debug)
    {
1096 1097 1098 1099
        Pout<< "globalMeshData::calcGlobalEdgeSlaves() :"
            << " coupled edges:" << edges.size()
            << " additional coupled edges:"
            << globalEdgeSlavesMapPtr_().constructSize() - edges.size()
mattijs's avatar
mattijs committed
1100 1101
            << endl;
    }
1102
}
mattijs's avatar
mattijs committed
1103 1104


1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122
void Foam::globalMeshData::calcGlobalEdgeOrientation() const
{
    if (debug)
    {
        Pout<< "globalMeshData::calcGlobalEdgeOrientation() :"
            << " calculating edge orientation w.r.t. master edge." << endl;
    }

    const globalIndex& globalPoints = globalPointNumbering();

    // 1. Determine master point
    labelList masterPoint;
    {
        const mapDistribute& map = globalPointSlavesMap();

        masterPoint.setSize(map.constructSize());
        masterPoint = labelMax;

1123
        for (label pointi = 0; pointi < coupledPatch().nPoints(); pointi++)
1124
        {
1125
            masterPoint[pointi] = globalPoints.toGlobal(pointi);
1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182
        }
        syncData
        (
            masterPoint,
            global