globalMeshData.C 77.7 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
6
     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
7
8
9
10
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

11
12
13
14
    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.
15
16
17
18
19
20
21

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

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

#include "globalMeshData.H"
#include "Time.H"
#include "Pstream.H"
#include "PstreamCombineReduceOps.H"
#include "processorPolyPatch.H"
#include "demandDrivenData.H"
#include "globalPoints.H"
33
34
#include "polyMesh.H"
#include "mapDistribute.H"
35
36
37
38
39
#include "labelIOList.H"
#include "PackedList.H"
#include "mergePoints.H"
#include "matchPoints.H"
#include "OFstream.H"
40
#include "globalIndexAndTransform.H"
41
42
43

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

44
45
46
namespace Foam
{
defineTypeNameAndDebug(globalMeshData, 0);
47
48

// Geometric matching tolerance. Factor of mesh bounding box.
49
const scalar globalMeshData::matchTol_ = 1e-8;
50

51
52
53
54
55
56
57
58
59
60
61
62
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]);
    }
};
}

63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

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

// Collect processor patch addressing.
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;

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


    if (Pstream::parRun())
    {
94
95
        PstreamBuffers pBufs(Pstream::nonBlocking);

96
        // Send indices of my processor patches to my neighbours
97
        forAll(processorPatches_, i)
98
99
100
        {
            label patchi = processorPatches_[i];

101
            UOPstream toNeighbour
102
103
104
105
            (
                refCast<const processorPolyPatch>
                (
                    mesh_.boundaryMesh()[patchi]
106
107
                ).neighbProcNo(),
                pBufs
108
109
110
111
112
            );

            toNeighbour << processorPatchIndices_[patchi];
        }

113
114
        pBufs.finishedSends();

115
116
117
        forAll(processorPatches_, i)
        {
            label patchi = processorPatches_[i];
118

119
            UIPstream fromNeighbour
120
121
122
123
            (
                refCast<const processorPolyPatch>
                (
                    mesh_.boundaryMesh()[patchi]
124
125
                ).neighbProcNo(),
                pBufs
126
            );
127

128
129
130
131
132
133
            fromNeighbour >> processorPatchNeighbours_[patchi];
        }
    }
}


134
135
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
246
247
248
249
250
251
252
253
void Foam::globalMeshData::calcSharedPoints() const
{
    if
    (
        nGlobalPoints_ != -1
     || sharedPointLabelsPtr_.valid()
     || sharedPointAddrPtr_.valid()
    )
    {
        FatalErrorIn("globalMeshData::calcSharedPoints()")
            << "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;
    }
}


254
255
256
// Given information about locally used edges allocate global shared edges.
void Foam::globalMeshData::countSharedEdges
(
257
258
    const EdgeMap<labelList>& procSharedEdges,
    EdgeMap<label>& globalShared,
259
260
261
262
    label& sharedEdgeI
)
{
    // Count occurrences of procSharedEdges in global shared edges table.
263
    forAllConstIter(EdgeMap<labelList>, procSharedEdges, iter)
264
265
266
    {
        const edge& e = iter.key();

267
        EdgeMap<label>::iterator globalFnd = globalShared.find(e);
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302

        if (globalFnd == globalShared.end())
        {
            // 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++);
            }
        }
        else
        {
            if (globalFnd() == -1)
            {
                // Second time occurence of this edge. Assign proper
                // edge label.
                globalFnd() = sharedEdgeI++;
            }
        }
    }
}


// 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)
void Foam::globalMeshData::calcSharedEdges() const
{
303
304
305
306
307
308
    if
    (
        nGlobalEdges_ != -1
     || sharedEdgeLabelsPtr_.valid()
     || sharedEdgeAddrPtr_.valid()
    )
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
    {
        FatalErrorIn("globalMeshData::calcSharedEdges()")
            << "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)
329
    EdgeMap<labelList> localShared(2*sharedPtAddr.size());
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354

    const edgeList& edges = mesh_.edges();

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

        Map<label>::const_iterator e0Fnd = meshToShared.find(e[0]);

        if (e0Fnd != meshToShared.end())
        {
            Map<label>::const_iterator e1Fnd = meshToShared.find(e[1]);

            if (e1Fnd != meshToShared.end())
            {
                // 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()]
                );

355
                EdgeMap<labelList>::iterator iter =
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
                    localShared.find(sharedEdge);

                if (iter == localShared.end())
                {
                    // 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();

                    label sz = edgeLabels.size();
                    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.

386
    EdgeMap<label> globalShared(nGlobalPoints());
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411

    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.
                IPstream fromSlave(Pstream::blocking, slave);
412
                EdgeMap<labelList> procSharedEdges(fromSlave);
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428

                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.
        {
429
            EdgeMap<label> oldSharedEdges(globalShared);
430
431
432

            globalShared.clear();

433
            forAllConstIter(EdgeMap<label>, oldSharedEdges, iter)
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
            {
                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.
                OPstream toSlave(Pstream::blocking, slave);
                toSlave << globalShared;
            }
        }
    }
    else
    {
        // Send local edges to master
        {
            OPstream toMaster(Pstream::blocking, Pstream::masterNo());

            toMaster << localShared;
        }
        // Receive merged edges from master.
        {
            IPstream fromMaster(Pstream::blocking, Pstream::masterNo());

            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());

489
    forAllConstIter(EdgeMap<labelList>, localShared, iter)
490
491
492
    {
        const edge& e = iter.key();

493
        EdgeMap<label>::const_iterator edgeFnd = globalShared.find(e);
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510

        if (edgeFnd != globalShared.end())
        {
            // My local edge is indeed a shared one. Go through all local edge
            // labels with this point combination.
            const labelList& edgeLabels = iter();

            forAll(edgeLabels, i)
            {
                // Store label of local mesh edge
                dynSharedEdgeLabels.append(edgeLabels[i]);

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

512
513
    sharedEdgeLabelsPtr_.reset(new labelList());
    labelList& sharedEdgeLabels = sharedEdgeLabelsPtr_();
514
515
    sharedEdgeLabels.transfer(dynSharedEdgeLabels);

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

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


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

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

543
    globalPointSlavesPtr_.reset
544
    (
545
546
547
548
        new labelListList
        (
            globalData.pointPoints().xfer()
        )
549
    );
550
    globalPointTransformedSlavesPtr_.reset
551
    (
552
        new labelListList
553
        (
554
            globalData.transformedPointPoints().xfer()
555
        )
556
    );
557

558
    globalPointSlavesMapPtr_.reset
559
560
561
    (
        new mapDistribute
        (
562
            globalData.map().xfer()
563
564
565
566
567
        )
    );
}


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


    // Create field with my local data
    labelPairList myData(globalPointSlavesMap().constructSize());
    forAll(slaves, pointI)
581
    {
582
583
584
585
586
587
        myData[pointI] = globalIndexAndTransform::encode
        (
            Pstream::myProcNo(),
            pointI,
            transforms.nullTransformIndex()
        );
588
    }
589
590
    // Send over.
    globalPointSlavesMap().distribute(myData);
591
592


593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
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
    // String of connected points with their transform
    allPointConnectivity.setSize(globalPointSlavesMap().constructSize());
    forAll(slaves, pointI)
    {
        // Reconstruct string of connected points
        const labelList& pSlaves = slaves[pointI];
        const labelList& pTransformSlaves = transformedSlaves[pointI];
        labelPairList& pConnectivity = allPointConnectivity[pointI];
        pConnectivity.setSize(1+pSlaves.size()+pTransformSlaves.size());
        label connI = 0;

        // Add myself
        pConnectivity[connI++] = myData[pointI];
        // 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]];
            label procI = globalIndexAndTransform::processor(n);
            label index = globalIndexAndTransform::index(n);
            pConnectivity[connI++] = globalIndexAndTransform::encode
            (
                procI,
                index,
                transformI
            );
        }

        // Put back in slots
        forAll(pSlaves, i)
        {
            allPointConnectivity[pSlaves[i]] = pConnectivity;
        }
        forAll(pTransformSlaves, i)
        {
            allPointConnectivity[pTransformSlaves[i]] = pConnectivity;
        }
    }
    globalPointSlavesMap().reverseDistribute
    (
        allPointConnectivity.size(),
        allPointConnectivity
    );
}
647
648


649
650
651
652
653
654
655
void Foam::globalMeshData::calcGlobalPointEdges
(
    labelListList& globalPointEdges,
    List<labelPairList>& globalPointPoints
) const
{
    const edgeList& edges = coupledPatch().edges();
656
    const labelListList& pointEdges = coupledPatch().pointEdges();
657
    const globalIndex& globalEdgeNumbers = globalEdgeNumbering();
658
659
660
661
662
663
    const labelListList& slaves = globalPointSlaves();
    const labelListList& transformedSlaves = globalPointTransformedSlaves();

    // Create local version
    globalPointEdges.setSize(globalPointSlavesMap().constructSize());
    globalPointPoints.setSize(globalPointSlavesMap().constructSize());
664
665
666
667
668
669
670
    forAll(pointEdges, pointI)
    {
        const labelList& pEdges = pointEdges[pointI];
        labelList& globalPEdges = globalPointEdges[pointI];
        globalPEdges.setSize(pEdges.size());
        forAll(pEdges, i)
        {
671
            globalPEdges[i] = globalEdgeNumbers.toGlobal(pEdges[i]);
672
        }
673
674
675
676
677
678
679
680
681
682
683
684
685

        labelPairList& globalPPoints = globalPointPoints[pointI];
        globalPPoints.setSize(pEdges.size());
        forAll(pEdges, i)
        {
            label otherPointI = edges[pEdges[i]].otherVertex(pointI);
            globalPPoints[i] = globalIndexAndTransform::encode
            (
                Pstream::myProcNo(),
                otherPointI,
                globalTransforms().nullTransformIndex()
            );
        }
686
687
    }

688
689
    // Pull slave data to master. Dummy transform.
    globalPointSlavesMap().distribute(globalPointEdges);
690
691
692
693
694
695
    globalPointSlavesMap().distribute(globalPointPoints);
    // Add all pointEdges
    forAll(slaves, pointI)
    {
        const labelList& pSlaves = slaves[pointI];
        const labelList& pTransformSlaves = transformedSlaves[pointI];
696

697
698
699
700
701
702
703
704
705
        label n = 0;
        forAll(pSlaves, i)
        {
            n += globalPointEdges[pSlaves[i]].size();
        }
        forAll(pTransformSlaves, i)
        {
            n += globalPointEdges[pTransformSlaves[i]].size();
        }
706

707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
        // Add all the point edges of the slaves to those of the (master) point
        {
            labelList& globalPEdges = globalPointEdges[pointI];
            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];
                }
            }
729

730
731
732
733
734
735
736
737
738
739
            // Put back in slots
            forAll(pSlaves, i)
            {
                globalPointEdges[pSlaves[i]] = globalPEdges;
            }
            forAll(pTransformSlaves, i)
            {
                globalPointEdges[pTransformSlaves[i]] = globalPEdges;
            }
        }
740

741

742
        // Same for corresponding pointPoints
743
        {
744
745
746
            labelPairList& globalPPoints = globalPointPoints[pointI];
            label sz = globalPPoints.size();
            globalPPoints.setSize(sz + n);
747

748
749
            // Add untransformed points
            forAll(pSlaves, i)
750
            {
751
752
                const labelPairList& otherData = globalPointPoints[pSlaves[i]];
                forAll(otherData, j)
753
                {
754
                    globalPPoints[sz++] = otherData[j];
755
756
                }
            }
757
758
            // Add transformed points.
            forAll(pTransformSlaves, i)
759
            {
760
761
762
763
764
                // Get transform from index
                label transformI = globalPointSlavesMap().whichTransform
                (
                    pTransformSlaves[i]
                );
765

766
767
768
769
770
771
772
773
774
775
776
777
778
779
                const labelPairList& otherData =
                    globalPointPoints[pTransformSlaves[i]];
                forAll(otherData, j)
                {
                    // Add transform to connectivity
                    const labelPair& n = otherData[j];
                    label procI = globalIndexAndTransform::processor(n);
                    label index = globalIndexAndTransform::index(n);
                    globalPPoints[sz++] = globalIndexAndTransform::encode
                    (
                        procI,
                        index,
                        transformI
                    );
780
781
                }
            }
782
783
784
785
786
787
788
789
790
791

            // Put back in slots
            forAll(pSlaves, i)
            {
                globalPointPoints[pSlaves[i]] = globalPPoints;
            }
            forAll(pTransformSlaves, i)
            {
                globalPointPoints[pTransformSlaves[i]] = globalPPoints;
            }
792
        }
793
794
795
796
797
798
799
800
801
802
803
804
805
806
    }
    // Push back
    globalPointSlavesMap().reverseDistribute
    (
        globalPointEdges.size(),
        globalPointEdges
    );
    // Push back
    globalPointSlavesMap().reverseDistribute
    (
        globalPointPoints.size(),
        globalPointPoints
    );
}
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
833
834
835
836
// Find transformation to take remotePoint to localPoint. Use info to find
// the transforms.
Foam::label Foam::globalMeshData::findTransform
(
    const labelPairList& info,
    const labelPair& remotePoint,
    const label localPoint
) const
{
    const label remoteProcI = globalIndexAndTransform::processor(remotePoint);
    const label remoteIndex = globalIndexAndTransform::index(remotePoint);

    label remoteTransformI = -1;
    label localTransformI = -1;
    forAll(info, i)
    {
        label procI = globalIndexAndTransform::processor(info[i]);
        label pointI = globalIndexAndTransform::index(info[i]);
        label transformI = globalIndexAndTransform::transformIndex(info[i]);

        if (procI == Pstream::myProcNo() && pointI == localPoint)
        {
            localTransformI = transformI;
            //Pout<< "For local :" << localPoint
            //    << " found transform:" << localTransformI
            //    << endl;
        }
        if (procI == remoteProcI && pointI == remoteIndex)
837
        {
838
839
840
841
842
843
844
            remoteTransformI = transformI;
            //Pout<< "For remote:" << remotePoint
            //    << " found transform:" << remoteTransformI
            //    << " at index:" << i
            //    << endl;
        }
    }
845

846
847
848
849
    if (remoteTransformI == -1 || localTransformI == -1)
    {
        FatalErrorIn("globalMeshData::findTransform(..)")
            << "Problem. Cannot find " << remotePoint
850
851
852
            << " or " << localPoint  << " "
            << coupledPatch().localPoints()[localPoint]
            << " in " << info
853
854
855
856
857
            << endl
            << "remoteTransformI:" << remoteTransformI << endl
            << "localTransformI:" << localTransformI
            << abort(FatalError);
    }
858

859
    return globalTransforms().subtractTransformIndex
860
861
862
863
864
    (
        remoteTransformI,
        localTransformI
    );
}
865

866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885

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();


    // The whole problem with deducting edge-connectivity from
    // point-connectivity is that one of the the endpoints might be
    // 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
886
    // 1. collect point connectivity - basically recreating globalPoints output.
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
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
    // All points will now have a string of points. The transforms are
    // 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
        (
            globalIndexAndTransform::encode
            (
                Pstream::myProcNo(),
                edgeI,
                globalTransforms().nullTransformIndex()
            )
        );

        forAll(pEdges0, i)
        {
            forAll(pEdges1, j)
936
            {
937
938
939
940
941
                if
                (
                    pEdges0[i] == pEdges1[j]
                 && pEdges0[i] != globalEdgeNumbers.toGlobal(edgeI)
                )
942
                {
943
944
945
946
947
                    // 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]

948

949
950
951
952
                    // Find difference in transforms to go from point on remote
                    // edge (pPoints1[j]) to this point.

                    label transform0 = findTransform
953
                    (
954
955
956
957
958
959
960
961
962
                        allPointConnectivity[e[0]],
                        pPoints1[j],
                        e[0]
                    );
                    label transform1 = findTransform
                    (
                        allPointConnectivity[e[1]],
                        pPoints0[i],
                        e[1]
963
964
                    );

965
966
967
968
                    if (transform0 == transform1)
                    {
                        label procI = globalEdgeNumbers.whichProcID(pEdges0[i]);
                        eEdges.append
969
970
971
972
                        (
                            globalIndexAndTransform::encode
                            (
                                procI,
973
974
                                globalEdgeNumbers.toLocal(procI, pEdges0[i]),
                                transform0
975
976
                            )
                        );
977
978
979
980
981
982
983
984
985
986
987
988
                    }
                }
            }
        }

        allEdgeConnectivity[edgeI].transfer(eEdges);
        sort(allEdgeConnectivity[edgeI], globalIndexAndTransform::less());
    }

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

990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
    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
            (
                (
                    globalIndexAndTransform::processor(masterInfo)
                 == Pstream::myProcNo()
                )
             && (globalIndexAndTransform::index(masterInfo) == edgeI)
            )
            {
                // 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];
                    label procI = globalIndexAndTransform::processor(info);
                    label index = globalIndexAndTransform::index(info);
                    label transform = globalIndexAndTransform::transformIndex
                    (
                        info
                    );

                    if (transform == globalTransforms().nullTransformIndex())
                    {
                        eEdges[nonTransformI++] = globalEdgeNumbers.toGlobal
                        (
                            procI,
                            index
                        );
                    }
                    else
                    {
                        trafoEEdges[transformI++] = info;
1041
1042
                    }
                }
1043
1044
1045

                eEdges.setSize(nonTransformI);
                trafoEEdges.setSize(transformI);
1046
1047
1048
1049
1050
            }
        }
    }


1051
1052
    // Construct map
    globalEdgeTransformedSlavesPtr_.reset(new labelListList());
1053
1054

    List<Map<label> > compactMap(Pstream::nProcs());
1055
    globalEdgeSlavesMapPtr_.reset
1056
1057
1058
    (
        new mapDistribute
        (
1059
1060
            globalEdgeNumbers,
            globalEdgeSlaves,
mattijs's avatar
mattijs committed
1061

1062
1063
1064
            globalTransforms(),
            transformedEdges,
            globalEdgeTransformedSlavesPtr_(),
mattijs's avatar
mattijs committed
1065

1066
1067
            compactMap
        )
mattijs's avatar
mattijs committed
1068
1069
1070
1071
1072
    );


    if (debug)
    {
1073
1074
1075
1076
        Pout<< "globalMeshData::calcGlobalEdgeSlaves() :"
            << " coupled edges:" << edges.size()
            << " additional coupled edges:"
            << globalEdgeSlavesMapPtr_().constructSize() - edges.size()
mattijs's avatar
mattijs committed
1077
1078
            << endl;
    }
1079
}
mattijs's avatar
mattijs committed
1080
1081


1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
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
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
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;

        for (label pointI = 0; pointI < coupledPatch().nPoints(); pointI++)
        {
            masterPoint[pointI] = globalPoints.toGlobal(pointI);
        }
        syncData
        (
            masterPoint,
            globalPointSlaves(),
            globalPointTransformedSlaves(),
            map,
            minEqOp<label>()
        );
    }

    // Now all points should know who is master by comparing their global
    // pointID with the masterPointID. We now can use this information
    // to find the orientation of the master edge.

    {
        const mapDistribute& map = globalEdgeSlavesMap();
        const labelListList& slaves = globalEdgeSlaves();
        const labelListList& transformedSlaves = globalEdgeTransformedSlaves();

        // Distribute orientation of master edge (in masterPoint numbering)
        labelPairList masterEdgeVerts(map.constructSize());
        masterEdgeVerts = labelPair(labelMax, labelMax);

        for (label edgeI = 0; edgeI < coupledPatch().nEdges(); edgeI++)
        {
            if
            (
                (
                    slaves[edgeI].size()
                  + transformedSlaves[edgeI].size()
                )
              > 0
            )
            {
                // I am master. Fill in my masterPoint equivalent.

                const edge& e = coupledPatch().edges()[edgeI];
                masterEdgeVerts[edgeI] = labelPair
                (
                    masterPoint[e[0]],
                    masterPoint[e[1]]
                );
            }
        }
        syncData
        (
            masterEdgeVerts,
            slaves,
            transformedSlaves,
            map,
            minEqOp<labelPair>()
        );

        // Now check my edges on how they relate to the master's edgeVerts
        globalEdgeOrientationPtr_.reset
        (
            new PackedBoolList(coupledPatch().nEdges())
        );
        PackedBoolList& globalEdgeOrientation = globalEdgeOrientationPtr_();

        forAll(coupledPatch().edges(), edgeI)
        {
            const edge& e = coupledPatch().edges()[edgeI];
            const labelPair masterE
            (
                masterPoint[e[0]],
                masterPoint[e[1]]
            );

            label stat = labelPair::compare
            (
                masterE,
                masterEdgeVerts[edgeI]
            );
            if (stat == 0)
            {
                FatalErrorIn
                (
                    "globalMeshData::calcGlobalEdgeOrientation() const"
                )   << "problem : my edge:" << e
                    << " in master points:" << masterE
                    << " v.s. masterEdgeVerts:" << masterEdgeVerts[edgeI]
                    << exit(FatalError);
            }
            else
            {
                globalEdgeOrientation[edgeI] = (stat == 1);
            }
        }
    }

    if (debug)
    {
        Pout<< "globalMeshData::calcGlobalEdgeOrientation() :"
            << " finished calculating edge orientation."
            << endl;
    }
}


1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
// Calculate uncoupled boundary faces (without calculating
// primitiveMesh::pointFaces())
void Foam::globalMeshData::calcPointBoundaryFaces
(
    labelListList& pointBoundaryFaces
) const
{
    const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
    const Map<label>& meshPointMap = coupledPatch().meshPointMap();

    // 1. Count

    labelList nPointFaces(coupledPatch().nPoints(), 0);

    forAll(bMesh, patchI)
    {
        const polyPatch& pp = bMesh[patchI];

        if (!pp.coupled())
        {
            forAll(pp, i)
            {
                const face& f = pp[i];

                forAll(f, fp)
                {
                    Map<label>::const_iterator iter = meshPointMap.find
                    (
                        f[fp]
                    );
                    if (iter != meshPointMap.end())
                    {
                        nPointFaces[iter()]++;
                    }
                }
            }
        }
    }


    // 2. Size

    pointBoundaryFaces.setSize(coupledPatch().nPoints());
    forAll(nPointFaces, pointI)
    {
        pointBoundaryFaces[pointI].setSize(nPointFaces[pointI]);
    }
    nPointFaces = 0;


    // 3. Fill

    forAll(bMesh, patchI)
    {
        const polyPatch& pp = bMesh[patchI];

        if (!pp.coupled())
        {
            forAll(pp, i)
            {
                const face& f = pp[i];
                forAll(f, fp)
                {
                    Map<label>::const_iterator iter = meshPointMap.find
                    (
                        f[fp]
                    );
                    if (iter != meshPointMap.end())
                    {
                        label bFaceI =
                             pp.start() + i - mesh_.nInternalFaces();
                        pointBoundaryFaces[iter()][nPointFaces[iter()]++] =
                            bFaceI;
                    }
                }
            }
        }
    }
}


void Foam::globalMeshData::calcGlobalPointBoundaryFaces() const
{
    if (debug)
    {
        Pout<< "globalMeshData::calcGlobalPointBoundaryFaces() :"
            << " calculating coupled point to boundary face addressing."
            << endl;
    }

    // Construct local point to (uncoupled)boundaryfaces.
    labelListList pointBoundaryFaces;
    calcPointBoundaryFaces(pointBoundaryFaces);


    // Global indices for boundary faces
    globalBoundaryFaceNumberingPtr_.reset
    (
        new globalIndex(mesh_.nFaces()-mesh_.nInternalFaces())
    );
    globalIndex& globalIndices = globalBoundaryFaceNumberingPtr_();


    // Convert local boundary faces to global numbering
    globalPointBoundaryFacesPtr_.reset
    (
        new labelListList(globalPointSlavesMap().constructSize())
    );
    labelListList& globalPointBoundaryFaces = globalPointBoundaryFacesPtr_();

    forAll(pointBoundaryFaces, pointI)
    {
        const labelList& bFaces = pointBoundaryFaces[pointI];
        labelList& globalFaces = globalPointBoundaryFaces[pointI];
        globalFaces.setSize(bFaces.size());
        forAll(bFaces, i)
        {
            globalFaces[i] = globalIndices.toGlobal(bFaces[i]);
        }
    }


    // Pull slave pointBoundaryFaces to master
    globalPointSlavesMap().distribute
    (
        globalPointBoundaryFaces,
        true    // put data on transformed points into correct slots
    );


    // Merge slave labels into master globalPointBoundaryFaces.
    // Split into untransformed and transformed values.
    const labelListList& pointSlaves = globalPointSlaves();
    const labelListList& pointTransformSlaves =
        globalPointTransformedSlaves();


    // Any faces coming in through transformation
    List<labelPairList> transformedFaces(pointSlaves.size());


    forAll(pointSlaves, pointI)
    {
        const labelList& slaves = pointSlaves[pointI];
        const labelList& transformedSlaves = pointTransformSlaves[pointI];

        if (slaves.size() > 0)
        {
            labelList& myBFaces = globalPointBoundaryFaces[pointI];
            label sz = myBFaces.size();

            // Count
            label n = 0;
            forAll(slaves, i)
            {
                n += globalPointBoundaryFaces[slaves[i]].size();
            }
            // Fill
            myBFaces.setSize(sz+n);
            n = sz;
            forAll(slaves, i)
            {
                const labelList& slaveBFaces =
                    globalPointBoundaryFaces[slaves[i]];

                // Add all slaveBFaces. Note that need to check for
                // uniqueness only in case of cyclics.

                forAll(slaveBFaces, j)
                {
                    label slave = slaveBFaces[j];
                    if (findIndex(SubList<label>(myBFaces, sz), slave) == -1)
                    {
                        myBFaces[n++] = slave;
                    }
                }
            }
            myBFaces.setSize(n);
        }


        if (transformedSlaves.size() > 0)
        {
            const labelList& untrafoFaces = globalPointBoundaryFaces[pointI];

            labelPairList& myBFaces = transformedFaces[pointI];
            label sz = myBFaces.size();

            // Count
            label n = 0;
            forAll(transformedSlaves, i)
            {
                n += globalPointBoundaryFaces[transformedSlaves[i]].size();
            }
            // Fill
            myBFaces.setSize(sz+n);
            n = sz;
            forAll(transformedSlaves, i)
            {
                label transformI = globalPointSlavesMap().whichTransform
                (
                    transformedSlaves[i]
                );

                const labelList& slaveBFaces =
                    globalPointBoundaryFaces[transformedSlaves[i]];

                forAll(slaveBFaces, j)
                {
                    label slave = slaveBFaces[j];
                    // Check that same face not already present untransformed
                    if (findIndex(untrafoFaces, slave)== -1)
                    {
                        label procI = globalIndices.whichProcID(slave);
                        label faceI = globalIndices.toLocal(procI, slave);

                        myBFaces[n++] = globalIndexAndTransform::encode
                        (
                            procI,
                            faceI,
                            transformI
                        );
                    }
                }
            }
            myBFaces.setSize(n);
        }


        if (slaves.size() + transformedSlaves.size() == 0)
        {
             globalPointBoundaryFaces[pointI].clear();
        }
    }

    // Construct a map to get the face data directly
    List<Map<label> > compactMap(Pstream::nProcs());

    globalPointTransformedBoundaryFacesPtr_.reset
    (
        new labelListList(transformedFaces.size())
    );

    globalPointBoundaryFacesMapPtr_.reset
    (
        new mapDistribute
        (
            globalIndices,
            globalPointBoundaryFaces,

            globalTransforms(),
            transformedFaces,
            globalPointTransformedBoundaryFacesPtr_(),

            compactMap
        )
    );
    globalPointBoundaryFaces.setSize(coupledPatch().nPoints());
    globalPointTransformedBoundaryFacesPtr_().setSize(coupledPatch().nPoints());

    if (debug)
    {
        Pout<< "globalMeshData::calcGlobalPointBoundaryFaces() :"
            << " coupled points:" << coupledPatch().nPoints()
            << " local boundary faces:" <<  globalIndices.localSize()
            << " additional coupled faces:"
            <<  globalPointBoundaryFacesMapPtr_().constructSize()
              - globalIndices.localSize()
            << endl;
    }
}


void Foam::globalMeshData::calcGlobalPointBoundaryCells() const
{
    if (debug)
    {
        Pout<< "globalMeshData::calcGlobalPointBoundaryCells() :"
            << " calculating coupled point to boundary cell addressing."
            << endl;
    }

    // Create map of boundary cells and point-cell addressing
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    label bCellI = 0;
    Map<label> meshCellMap(4*coupledPatch().nPoints());
    DynamicList<label> cellMap(meshCellMap.size());

    // Create addressing for point to boundary cells (local)
    labelListList pointBoundaryCells(coupledPatch().nPoints());

    forAll(coupledPatch().meshPoints(), pointI)
    {
        label meshPointI = coupledPatch().meshPoints()[pointI];
        const labelList& pCells = mesh_.pointCells(meshPointI);

        labelList& bCells = pointBoundaryCells[pointI];
        bCells.setSize(pCells.size());

        forAll(pCells, i)
        {
            label cellI = pCells[i];
            Map<label>::iterator fnd = meshCellMap.find(cellI);

            if (fnd != meshCellMap.end())
            {
                bCells[i] = fnd();
            }
            else
            {
                meshCellMap.insert(cellI, bCellI);
                cellMap.append(cellI);
                bCells[i] = bCellI;
                bCellI++;
            }
        }
    }


    boundaryCellsPtr_.reset(new labelList());
    labelList& boundaryCells = boundaryCellsPtr_();
    boundaryCells.transfer(cellMap.shrink());


    // Convert point-cells to global (boundary)cell numbers
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    globalBoundaryCellNumberingPtr_.reset
    (
        new globalIndex(boundaryCells.size())
    );
    globalIndex& globalIndices = globalBoundaryCellNumberingPtr_();


    globalPointBoundaryCellsPtr_.reset
    (
        new labelListList(globalPointSlavesMap().constructSize())
    );
    labelListList& globalPointBoundaryCells = globalPointBoundaryCellsPtr_();

    forAll(pointBoundaryCells, pointI)
    {
        const labelList& pCells = pointBoundaryCells[pointI];
        labelList& globalCells = globalPointBoundaryCells[pointI];
        globalCells.setSize(pCells.size());
        forAll(pCells, i)
        {
            globalCells[i] = globalIndices.toGlobal(pCells[i]);
        }
    }


    // Pull slave pointBoundaryCells to master
    globalPointSlavesMap().distribute
    (
        globalPointBoundaryCells,
        true    // put data on transformed points into correct slots
    );


    // Merge slave labels into master globalPointBoundaryCells
    const labelListList& pointSlaves = globalPointSlaves();
    const labelListList& pointTransformSlaves =
        globalPointTransformedSlaves();

    List<labelPairList> transformedCells(pointSlaves.size());


    forAll(pointSlaves, pointI)
    {
        const labelList& slaves = pointSlaves[pointI];
        const labelList& transformedSlaves = pointTransformSlaves[pointI];

        if (slaves.size() > 0)
        {
            labelList& myBCells = globalPointBoundaryCells[pointI];
            label sz = myBCells.size();

            // Count
            label n = 0;
            forAll(slaves, i)
            {
                n += globalPointBoundaryCells[slaves[i]].size();
            }
            // Fill
            myBCells.setSize(sz+n);
            n = sz;
            forAll(slaves, i)
            {
                const labelList& slaveBCells =
                    globalPointBoundaryCells[slaves[i]];

                // Add all slaveBCells. Note that need to check for
                // uniqueness only in case of cyclics.

                forAll(slaveBCells, j)
                {
                    label slave = slaveBCells[j];
                    if (findIndex(SubList<label>(myBCells, sz), slave) == -1)
                    {
                        myBCells[n++] = slave;
                    }
                }
            }
            myBCells.setSize(n);
        }


        if (transformedSlaves.size() > 0)
        {
            const labelList& untrafoCells = globalPointBoundaryCells[pointI];

            labelPairList& myBCells = transformedCells[pointI];
            label sz = myBCells.size();

            // Count
            label n = 0;
            forAll(transformedSlaves, i)
            {
                n += globalPointBoundaryCells[transformedSlaves[i]].size();
            }
            // Fill
            myBCells.setSize(sz+n);
            n = sz;
            forAll(transformedSlaves, i)
            {
                label transformI = globalPointSlavesMap().whichTransform
                (
                    transformedSlaves[i]
                );

                const labelList& slaveBCells =
                    globalPointBoundaryCells[transformedSlaves[i]];

                forAll(slaveBCells, j)
                {
                    label slave = slaveBCells[j];

                    // Check that same cell not already present untransformed
                    if (findIndex(untrafoCells, slave)== -1)
                    {
                        label procI = globalIndices.whichProcID(slave);
                        label cellI = globalIndices.toLocal(procI, slave);
                        myBCells[n++] = globalIndexAndTransform::encode
                        (
                            procI,
                            cellI,
                            transformI
                        );
                    }
                }
            }
            myBCells.setSize(n);
        }

        if (slaves.size() + transformedSlaves.size() == 0)
        {
             globalPointBoundaryCells[pointI].clear();
        }
    }

    // Construct a map to get the cell data directly
    List<Map<label> > compactMap(Pstream::nProcs());

    globalPointTransformedBoundaryCellsPtr_.reset
    (
        new labelListList(transformedCells.size())
    );

    globalPointBoundaryCellsMapPtr_.reset
    (
        new mapDistribute
        (
            globalIndices,
            globalPointBoundaryCells,

            globalTransforms(),
            transformedCells,
            globalPointTransformedBoundaryCellsPtr_(),

            compactMap
        )
    );
    globalPointBoundaryCells.setSize(coupledPatch().nPoints());
    globalPointTransformedBoundaryCellsPtr_().setSize(coupledPatch().nPoints());

    if (debug)
    {
        Pout<< "globalMeshData::calcGlobalPointBoundaryCells() :"
            << " coupled points:" << coupledPatch().nPoints()
            << " local boundary cells:" <<  globalIndices.localSize()
            << " additional coupled cells:"
            <<  globalPointBoundaryCellsMapPtr_().constructSize()
              - globalIndices.localSize()
            << endl;
    }
}
mattijs's avatar
mattijs committed
1702
1703


1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
void Foam::globalMeshData::calcGlobalCoPointSlaves() const
{
    if (debug)
    {
        Pout<< "globalMeshData::calcGlobalCoPointSlaves() :"
            << " calculating coupled master to collocated"
            << " slave point addressing." << endl;
    }

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

    globalCoPointSlavesPtr_.reset
    (
        new labelListList
        (
            globalData.pointPoints().xfer()
        )
    );
    globalCoPointSlavesMapPtr_.reset
    (
        new mapDistribute
        (
            globalData.map().xfer()
        )
    );

    if (debug)
    {
        Pout<< "globalMeshData::calcGlobalCoPointSlaves() :"
            << " finished calculating coupled master to collocated"
            << " slave point addressing." << endl;
    }
}


1740
1741
1742
1743
1744
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

// Construct from polyMesh
Foam::globalMeshData::globalMeshData(const polyMesh& mesh)
:
1745
    processorTopology(mesh.boundaryMesh(), UPstream::worldComm),
1746
1747
1748
1749
1750
1751
1752
1753
    mesh_(mesh),
    nTotalPoints_(-1),
    nTotalFaces_(-1),
    nTotalCells_(-1),
    processorPatches_(0),
    processorPatchIndices_(0),
    processorPatchNeighbours_(0),
    nGlobalPoints_(-1),
1754
1755
    sharedPointLabelsPtr_(NULL),
    sharedPointAddrPtr_(NULL),
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
    sharedPointGlobalLabelsPtr_(NULL),
    nGlobalEdges_(-1),
    sharedEdgeLabelsPtr_(NULL),
    sharedEdgeAddrPtr_(NULL)
{
    updateMesh();
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::globalMeshData::~globalMeshData()
{
    clearOut();
}


void Foam::globalMeshData::clearOut()
{
1775
    // Point
1776
    nGlobalPoints_ = -1;
1777
1778
1779
    sharedPointLabelsPtr_.clear();
    sharedPointAddrPtr_.clear();
    sharedPointGlobalLabelsPtr_.clear();
1780

1781
1782
1783
1784
1785
1786
    // Edge
    nGlobalEdges_ = -1;
    sharedEdgeLabelsPtr_.clear();
    sharedEdgeAddrPtr_.clear();

    // Coupled patch
1787
    coupledPatchPtr_.clear();
1788
1789
    coupledPatchMeshEdgesPtr_.clear();
    coupledPatchMeshEdgeMapPtr_.clear();
1790
    globalTransformsPtr_.clear();
1791

1792
1793
1794
    // Point
    globalPointNumberingPtr_.clear();
    globalPointSlavesPtr_.clear();
1795
    globalPointTransformedSlavesPtr_.clear();
1796
1797
1798
1799
    globalPointSlavesMapPtr_.clear();
    // Edge
    globalEdgeNumberingPtr_.clear();
    globalEdgeSlavesPtr_.clear();
1800
    globalEdgeTransformedSlavesPtr_.clear();
1801
    globalEdgeOrientationPtr_.clear();
1802
    globalEdgeSlavesMapPtr_.clear();
mattijs's avatar
mattijs committed
1803

1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
    // Face
    globalBoundaryFaceNumberingPtr_.clear();
    globalPointBoundaryFacesPtr_.clear();
    globalPointTransformedBoundaryFacesPtr_.clear();
    globalPointBoundaryFacesMapPtr_.clear();

    // Cell
    boundaryCellsPtr_.clear();
    globalBoundaryCellNumberingPtr_.clear();
    globalPointBoundaryCellsPtr_.clear();
    globalPointTransformedBoundaryCellsPtr_.clear();
    globalPointBoundaryCellsMapPtr_.clear();
1816
1817
1818
1819

    // Other: collocated points
    globalCoPointSlavesPtr_.clear();
    globalCoPointSlavesMapPtr_.clear();
1820
1821
1822
1823
1824
1825
1826
1827
}


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

// Return shared point global labels.
const Foam::labelList& Foam::globalMeshData::sharedPointGlobalLabels() const
{
1828
    if (!sharedPointGlobalLabelsPtr_.valid())
1829
    {
1830
1831
1832
1833
1834
        sharedPointGlobalLabelsPtr_.reset
        (
            new labelList(sharedPointLabels().size())
        );
        labelList& sharedPointGlobalLabels = sharedPointGlobalLabelsPtr_();
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852

        IOobject addrHeader
        (
            "pointProcAddressing",
            mesh_.facesInstance()/mesh_.meshSubDir,
            mesh_,
            IOobject::MUST_READ
        );

        if (addrHeader.headerOk())
        {
            // There is a pointProcAddressing file so use it to get labels
            // on the original mesh
            Pout<< "globalMeshData::sharedPointGlobalLabels : "
                << "Reading pointProcAddressing" << endl;

            labelIOList pointProcAddressing(addrHeader);

1853
1854
1855
            const labelList& pointLabels = sharedPointLabels();

            forAll(pointLabels, i)
1856
1857
            {
                // Get my mesh point
1858
                label pointI = pointLabels[i];
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871

                // Map to mesh point of original mesh
                sharedPointGlobalLabels[i] = pointProcAddressing[pointI];
            }
        }
        else
        {
            Pout<< "globalMeshData::sharedPointGlobalLabels :"
                << " Setting pointProcAddressing to -1" << endl;

            sharedPointGlobalLabels = -1;
        }
    }
1872
    return sharedPointGlobalLabelsPtr_();
1873
1874
1875
1876
1877
1878
1879
1880
1881
}


// Collect coordinates of shared points. (does parallel communication!)
Foam::pointField Foam::globalMeshData::sharedPoints() const
{
    // Get all processors to send their shared points to master.
    // (not very efficient)

1882
1883
1884
    pointField sharedPoints(nGlobalPoints());
    const labelList& pointAddr = sharedPointAddr();
    const labelList& pointLabels = sharedPointLabels();
1885
1886
1887
1888
1889

    if (Pstream::master())
    {
        // Master:
        // insert my own data first
1890
        forAll(pointLabels, i)
1891
        {
1892
            label sharedPointI = pointAddr[i];
1893

1894
            sharedPoints[sharedPointI] = mesh_.points()[pointLabels[i]];
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
        }

        // Receive data from slaves and insert
        for
        (
            int slave=Pstream::firstSlave();
            slave<=Pstream::lastSlave();
            slave++
        )
        {
            IPstream fromSlave(Pstream::blocking, slave);

            labelList nbrSharedPointAddr;
            pointField nbrSharedPoints;
            fromSlave >> nbrSharedPointAddr >> nbrSharedPoints;

            forAll(nbrSharedPointAddr, i)
            {
                label sharedPointI = nbrSharedPointAddr[i];

                sharedPoints[sharedPointI] = nbrSharedPoints[i];
            }
        }

        // Send back
        for
        (
            int slave=Pstream::firstSlave();
            slave<=Pstream::lastSlave();
            slave++
        )
        {
            OPstream toSlave
            (
                Pstream::blocking,
                slave,
                sharedPoints.size()*sizeof(vector::zero)
            );
            toSlave << sharedPoints;
        }
    }
    else
    {
        // Slave:
        // send points
        {
            OPstream toMaster(Pstream::blocking, Pstream::masterNo());

            toMaster
1944
1945
                << pointAddr
                << UIndirectList<point>(mesh_.points(), pointLabels)();
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
        }

        // Receive sharedPoints
        {
            IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
            fromMaster >> sharedPoints;
        }
    }

    return sharedPoints;
}


// Collect coordinates of shared points. (does parallel communication!)
Foam::pointField Foam::globalMeshData::geometricSharedPoints() const
{
    // Get coords of my shared points
1963
    pointField sharedPoints(mesh_.points(), sharedPointLabels());
1964
1965

    // Append from all processors
1966
    combineReduce(sharedPoints, ListPlusEqOp<pointField>());
1967
1968

    // Merge tolerance
1969
    scalar tolDim = matchTol_ * mesh_.bounds().mag();
1970
1971
1972
1973
1974

    // And see how many are unique
    labelList pMap;
    pointField mergedPoints;

mattijs's avatar
mattijs committed
1975
    Foam::mergePoints
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
    (
        sharedPoints,   // coordinates to merge
        tolDim,         // tolerance
        false,          // verbosity
        pMap,
        mergedPoints
    );

    return mergedPoints;
}


1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
Foam::label Foam::globalMeshData::nGlobalPoints() const
{
    if (nGlobalPoints_ == -1)
    {
        calcSharedPoints();
    }
    return nGlobalPoints_;
}


const Foam::labelList& Foam::globalMeshData::sharedPointLabels() const
{
    if (!sharedPointLabelsPtr_.valid())
    {
        calcSharedPoints();
    }
    return sharedPointLabelsPtr_();
}


const Foam::labelList& Foam::globalMeshData::sharedPointAddr() const
{
    if (!sharedPointAddrPtr_.valid())
    {
        calcSharedPoints();
    }
    return sharedPointAddrPtr_();
}


2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
Foam::label Foam::globalMeshData::nGlobalEdges() const
{
    if (nGlobalEdges_ == -1)
    {
        calcSharedEdges();
    }
    return nGlobalEdges_;
}


const Foam::labelList& Foam::globalMeshData::sharedEdgeLabels() const
{
2030
    if (!sharedEdgeLabelsPtr_.valid())
2031
2032
2033
    {
        calcSharedEdges();
    }
2034
    return sharedEdgeLabelsPtr_();
2035
2036
2037
2038
2039
}


const Foam::labelList& Foam::globalMeshData::sharedEdgeAddr() const
{
2040
    if (!sharedEdgeAddrPtr_.valid())