mappedPatchBase.C 54.5 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
OpenFOAM bot's avatar
OpenFOAM bot committed
6
7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2016 OpenFOAM Foundation
9
    Copyright (C) 2015-2020 OpenCFD Ltd.
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

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

#include "mappedPatchBase.H"
#include "addToRunTimeSelectionTable.H"
#include "ListListOps.H"
32
#include "meshSearchMeshObject.H"
33
34
35
36
#include "meshTools.H"
#include "OFstream.H"
#include "Random.H"
#include "treeDataFace.H"
37
#include "treeDataPoint.H"
38
39
40
41
42
43
#include "indexedOctree.H"
#include "polyMesh.H"
#include "polyPatch.H"
#include "Time.H"
#include "mapDistribute.H"
#include "SubField.H"
44
#include "triPointRef.H"
45
#include "syncTools.H"
46
#include "treeDataCell.H"
47
#include "DynamicField.H"
Andrew Heather's avatar
Andrew Heather committed
48
#include "faceAreaWeightAMI.H"
49
#include "OTstream.H"
50
51
52
53
54
55
56
57
58

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

namespace Foam
{
    defineTypeNameAndDebug(mappedPatchBase, 0);
}


59
60
61
62
63
const Foam::Enum
<
    Foam::mappedPatchBase::sampleMode
>
Foam::mappedPatchBase::sampleModeNames_
Mark OLESEN's avatar
Mark OLESEN committed
64
({
65
66
67
68
69
70
    { sampleMode::NEARESTCELL, "nearestCell" },
    { sampleMode::NEARESTPATCHFACE, "nearestPatchFace" },
    { sampleMode::NEARESTPATCHFACEAMI, "nearestPatchFaceAMI" },
    { sampleMode::NEARESTPATCHPOINT, "nearestPatchPoint" },
    { sampleMode::NEARESTFACE, "nearestFace" },
    { sampleMode::NEARESTONLYCELL, "nearestOnlyCell" },
Mark OLESEN's avatar
Mark OLESEN committed
71
});
72
73
74
75
76
77
78


const Foam::Enum
<
    Foam::mappedPatchBase::offsetMode
>
Foam::mappedPatchBase::offsetModeNames_
Mark OLESEN's avatar
Mark OLESEN committed
79
({
80
81
82
    { offsetMode::UNIFORM, "uniform" },
    { offsetMode::NONUNIFORM, "nonuniform" },
    { offsetMode::NORMAL, "normal" },
Mark OLESEN's avatar
Mark OLESEN committed
83
});
84
85
86
87


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

88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
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
Foam::autoPtr<Foam::fileName> Foam::mappedPatchBase::readDatabase
(
    const dictionary& dict
)
{
    autoPtr<fileName> dbNamePtr_;

    if (dict.found("sampleDatabase"))
    {
        const bool useDb = dict.get<bool>("sampleDatabase");
        if (useDb)
        {
            dbNamePtr_.set
            (
                new fileName
                (
                    dict.lookupOrDefault<fileName>
                    (
                        "sampleDatabasePath",
                        fileName::null
                    )
                )
            );
        }
    }
    else if (dict.found("sampleDatabasePath"))
    {
        dbNamePtr_.set(new fileName(dict.get<fileName>("sampleDatabasePath")));
    }

    return dbNamePtr_;
}


Foam::label Foam::mappedPatchBase::communicator
(
    const word& sampleWorld
)
{
    // Start off with local world
    label comm = UPstream::worldComm;

    if (!sampleWorld.empty() && Pstream::parRun())
    {
        if (!UPstream::allWorlds().found(sampleWorld))
        {
            FatalErrorInFunction << "Cannot find sampleWorld " << sampleWorld
                << " in set of worlds " << UPstream::allWorlds()
                << exit(FatalError);
        }

        const labelList& worldIDs = UPstream::worldIDs();

        DynamicList<label> subRanks(worldIDs.size());
        forAll(worldIDs, proci)
        {
            const label worldi = worldIDs[proci];
            if
            (
                worldi == UPstream::myWorldID()
             || UPstream::allWorlds()[worldi] == sampleWorld
            )
            {
                subRanks.append(proci);
            }
        }

        // Allocate new communicator with parent 0 (= world)
        comm = UPstream::allocateCommunicator(0, subRanks, true);

        if (debug)
        {
            Pout<< "mappedPatchBase::communicator :"
                << " myWorld:" << UPstream::myWorld()
                << " sampleWorld:" << sampleWorld
                << " using subRanks:" << subRanks
                << " new comm:" << comm << endl;
        }
    }

    return comm;
}


172
173
174
175
176
177
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::facePoints
(
    const polyPatch& pp
) const
{
    const polyMesh& mesh = pp.boundaryMesh().mesh();
178
179
180

    // Force construction of min-tet decomp
    (void)mesh.tetBasePtIs();
181
182

    // Initialise to face-centre
183
    tmp<pointField> tfacePoints(new pointField(patch_.size()));
184
    pointField& facePoints = tfacePoints.ref();
185

186
    forAll(pp, facei)
187
    {
188
        facePoints[facei] = facePoint
189
190
        (
            mesh,
191
            pp.start()+facei,
192
            polyMesh::FACE_DIAG_TRIS
193
        ).rawPoint();
194
195
196
197
198
199
    }

    return tfacePoints;
}


200
201
void Foam::mappedPatchBase::collectSamples
(
202
    const label mySampleWorld,      // Wanted world
203
    const pointField& facePoints,
204
205
    pointField& samples,            // All samples
    labelList& patchFaceWorlds,     // Per sample: wanted world
Mark OLESEN's avatar
Mark OLESEN committed
206
    labelList& patchFaceProcs,      // Per sample: originating processor
207
208
    labelList& patchFaces,          // Per sample: originating patchFace index
    pointField& patchFc             // Per sample: originating centre
209
210
) const
{
211
212
213
214
215
216
    const label oldComm(Pstream::warnComm);
    Pstream::warnComm = comm_;

    const label myRank = Pstream::myProcNo(comm_);
    const label nProcs = Pstream::nProcs(comm_);

217
    // Collect all sample points and the faces they come from.
218
    {
219
220
221
222
        List<pointField> globalFc(nProcs);
        globalFc[myRank] = facePoints;
        Pstream::gatherList(globalFc, Pstream::msgType(), comm_);
        Pstream::scatterList(globalFc, Pstream::msgType(), comm_);
223
224
225
226
227
228
229
        // Rework into straight list
        patchFc = ListListOps::combine<pointField>
        (
            globalFc,
            accessOp<pointField>()
        );
    }
230

231
    {
232
233
234
235
        List<pointField> globalSamples(nProcs);
        globalSamples[myRank] = samplePoints(facePoints);
        Pstream::gatherList(globalSamples, Pstream::msgType(), comm_);
        Pstream::scatterList(globalSamples, Pstream::msgType(), comm_);
236
237
238
239
240
241
242
243
244
        // Rework into straight list
        samples = ListListOps::combine<pointField>
        (
            globalSamples,
            accessOp<pointField>()
        );
    }

    {
245
246
        labelListList globalFaces(nProcs);
        globalFaces[myRank] = identity(patch_.size());
247
        // Distribute to all processors
248
249
        Pstream::gatherList(globalFaces, Pstream::msgType(), comm_);
        Pstream::scatterList(globalFaces, Pstream::msgType(), comm_);
250
251

        patchFaces = ListListOps::combine<labelList>
252
253
254
        (
            globalFaces,
            accessOp<labelList>()
255
256
257
        );
    }

258
    {
259
260
261
262
263
264
265
266
267
        labelList procToWorldIndex(nProcs);
        procToWorldIndex[myRank] = mySampleWorld;
        Pstream::gatherList(procToWorldIndex, Pstream::msgType(), comm_);
        Pstream::scatterList(procToWorldIndex, Pstream::msgType(), comm_);

        labelList nPerProc(nProcs);
        nPerProc[myRank] = patch_.size();
        Pstream::gatherList(nPerProc, Pstream::msgType(), comm_);
        Pstream::scatterList(nPerProc, Pstream::msgType(), comm_);
268

269
        patchFaceWorlds.setSize(patchFaces.size());
270
271
272
        patchFaceProcs.setSize(patchFaces.size());

        label sampleI = 0;
273
        forAll(nPerProc, proci)
274
        {
275
            for (label i = 0; i < nPerProc[proci]; i++)
276
            {
277
278
279
                patchFaceWorlds[sampleI] = procToWorldIndex[proci];
                patchFaceProcs[sampleI] = proci;
                sampleI++;
280
            }
281
282
        }
    }
283
    Pstream::warnComm = oldComm;
284
285
286
}


287
void Foam::mappedPatchBase::findLocalSamples
288
(
289
    const sampleMode mode,
290
291
292
293
294

    const label mySampleWorld,  // local world to sample == my own world
    const word& sampleRegion,   // local region to sample
    const word& samplePatch,    // local patch to sample

295
    const pointField& samples,
296
    List<nearInfoWorld>& nearest
297
298
) const
{
299
300
301
302
    // Find the local cell containing the samples

    const label myRank = Pstream::myProcNo(comm_);

303
    // Lookup the correct region
304
    const polyMesh& mesh = lookupMesh(sampleRegion);
305
306

    // All the info for nearest. Construct to miss
307
308
309
310
311
312
313
    nearest.setSize(samples.size());
    nearInfoWorld miss;
    {
        miss.first().second() = Tuple2<scalar, label>(Foam::sqr(GREAT), -1);
        miss.second() = -1; // set world to be ignored
    }
    nearest = miss;
314

315
    switch (mode)
316
317
318
    {
        case NEARESTCELL:
        {
319
            if (samplePatch.size() && samplePatch != "none")
320
            {
321
322
                FatalErrorInFunction
                    << "No need to supply a patch name when in "
323
                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
324
325
            }

326
            //- Note: face-diagonal decomposition
327
            const indexedOctree<Foam::treeDataCell>& tree = mesh.cellTree();
328
329
330
331

            forAll(samples, sampleI)
            {
                const point& sample = samples[sampleI];
332
                nearInfoWorld& near = nearest[sampleI];
333

334
                label celli = tree.findInside(sample);
335

336
                if (celli == -1)
337
                {
338
339
340
                    near.first().second().first() = Foam::sqr(GREAT);
                    near.first().second().second() = myRank;
                    near.second() = mySampleWorld;
341
342
343
                }
                else
                {
344
                    const point& cc = mesh.cellCentres()[celli];
345

346
                    near.first().first() = pointIndexHit
347
348
349
                    (
                        true,
                        cc,
350
                        celli
351
                    );
352
353
354
                    near.first().second().first() = magSqr(cc-sample);
                    near.first().second().second() = myRank;
                    near.second() = mySampleWorld;
355
356
357
358
359
                }
            }
            break;
        }

360
361
        case NEARESTONLYCELL:
        {
362
            if (samplePatch.size() && samplePatch != "none")
363
            {
364
365
                FatalErrorInFunction
                    << "No need to supply a patch name when in "
366
367
368
369
370
371
372
373
374
                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
            }

            //- Note: face-diagonal decomposition
            const indexedOctree<Foam::treeDataCell>& tree = mesh.cellTree();

            forAll(samples, sampleI)
            {
                const point& sample = samples[sampleI];
375
                nearInfoWorld& near = nearest[sampleI];
376

377
378
                near.first().first() = tree.findNearest(sample, sqr(GREAT));
                near.first().second().first() = magSqr
379
                (
380
                    near.first().first().hitPoint()
381
382
                   -sample
                );
383
384
                near.first().second().second() = myRank;
                near.second() = mySampleWorld;
385
386
387
388
            }
            break;
        }

389
390
391
392
        case NEARESTPATCHFACE:
        {
            Random rndGen(123456);

393
            const polyPatch& pp = lookupPatch(sampleRegion, samplePatch);
394
395
396
397
398

            if (pp.empty())
            {
                forAll(samples, sampleI)
                {
399
400
401
402
                    nearInfoWorld& near = nearest[sampleI];
                    near.first().second().first() = Foam::sqr(GREAT);
                    near.first().second().second() = myRank;
                    near.second() = mySampleWorld;
403
404
405
406
407
408
409
410
411
                }
            }
            else
            {
                treeBoundBox patchBb
                (
                    treeBoundBox(pp.points(), pp.meshPoints()).extend
                    (
                        rndGen,
412
                        1e-4
413
414
                    )
                );
415
416
                patchBb.min() -= point::uniform(ROOTVSMALL);
                patchBb.max() += point::uniform(ROOTVSMALL);
417
418
419
420
421
422
423

                indexedOctree<treeDataFace> boundaryTree
                (
                    treeDataFace    // all information needed to search faces
                    (
                        false,      // do not cache bb
                        mesh,
424
                        identity(pp.range())  // boundary faces only
425
426
427
428
429
430
431
432
433
434
435
                    ),
                    patchBb,        // overall search domain
                    8,              // maxLevel
                    10,             // leafsize
                    3.0             // duplicity
                );

                forAll(samples, sampleI)
                {
                    const point& sample = samples[sampleI];

436
437
                    nearInfoWorld& near = nearest[sampleI];
                    pointIndexHit& nearInfo = near.first().first();
438
439
440
441
442
443
444
445
                    nearInfo = boundaryTree.findNearest
                    (
                        sample,
                        magSqr(patchBb.span())
                    );

                    if (!nearInfo.hit())
                    {
446
447
448
                        near.first().second().first() = Foam::sqr(GREAT);
                        near.first().second().second() = myRank;
                        near.second() = mySampleWorld;
449
450
451
452
453
                    }
                    else
                    {
                        point fc(pp[nearInfo.index()].centre(pp.points()));
                        nearInfo.setPoint(fc);
454
455
456
                        near.first().second().first() = magSqr(fc-sample);
                        near.first().second().second() = myRank;
                        near.second() = mySampleWorld;
457
458
459
460
461
462
                    }
                }
            }
            break;
        }

463
464
465
466
        case NEARESTPATCHPOINT:
        {
            Random rndGen(123456);

467
            const polyPatch& pp = lookupPatch(sampleRegion, samplePatch);
468
469
470
471
472

            if (pp.empty())
            {
                forAll(samples, sampleI)
                {
473
474
475
476
                    nearInfoWorld& near = nearest[sampleI];
                    near.first().second().first() = Foam::sqr(GREAT);
                    near.first().second().second() = myRank;
                    near.second() = mySampleWorld;
477
478
479
480
481
482
483
484
485
486
487
488
489
                }
            }
            else
            {
                // patch (local) points
                treeBoundBox patchBb
                (
                    treeBoundBox(pp.points(), pp.meshPoints()).extend
                    (
                        rndGen,
                        1e-4
                    )
                );
490
491
                patchBb.min() -= point::uniform(ROOTVSMALL);
                patchBb.max() += point::uniform(ROOTVSMALL);
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509

                indexedOctree<treeDataPoint> boundaryTree
                (
                    treeDataPoint   // all information needed to search faces
                    (
                        mesh.points(),
                        pp.meshPoints() // selection of points to search on
                    ),
                    patchBb,        // overall search domain
                    8,              // maxLevel
                    10,             // leafsize
                    3.0             // duplicity
                );

                forAll(samples, sampleI)
                {
                    const point& sample = samples[sampleI];

510
511
                    nearInfoWorld& near = nearest[sampleI];
                    pointIndexHit& nearInfo = near.first().first();
512
513
514
515
516
517
518
519
                    nearInfo = boundaryTree.findNearest
                    (
                        sample,
                        magSqr(patchBb.span())
                    );

                    if (!nearInfo.hit())
                    {
520
521
522
                        near.first().second().first() = Foam::sqr(GREAT);
                        near.first().second().second() = myRank;
                        near.second() = mySampleWorld;
523
524
525
526
527
                    }
                    else
                    {
                        const point& pt = nearInfo.hitPoint();

528
529
530
                        near.first().second().first() = magSqr(pt-sample);
                        near.first().second().second() = myRank;
                        near.second() = mySampleWorld;
531
532
533
534
535
536
                    }
                }
            }
            break;
        }

537
538
        case NEARESTFACE:
        {
539
            if (samplePatch.size() && samplePatch != "none")
540
            {
541
542
                FatalErrorInFunction
                    << "No need to supply a patch name when in "
543
                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
544
545
            }

546
547
548
            //- Note: face-diagonal decomposition
            const meshSearchMeshObject& meshSearchEngine =
                meshSearchMeshObject::New(mesh);
549
550
551
552

            forAll(samples, sampleI)
            {
                const point& sample = samples[sampleI];
553
                nearInfoWorld& near = nearest[sampleI];
554

555
                label facei = meshSearchEngine.findNearestFace(sample);
556

557
                if (facei == -1)
558
                {
559
560
561
                    near.first().second().first() = Foam::sqr(GREAT);
                    near.first().second().second() = myRank;
                    near.second() = mySampleWorld;
562
563
564
                }
                else
                {
565
                    const point& fc = mesh.faceCentres()[facei];
566

567
568
569
570
                    near.first().first() = pointIndexHit(true, fc, facei);
                    near.first().second().first() = magSqr(fc-sample);
                    near.first().second().second() = myRank;
                    near.second() = mySampleWorld;
571
572
573
574
575
576
577
578
579
580
581
582
583
                }
            }
            break;
        }

        case NEARESTPATCHFACEAMI:
        {
            // nothing to do here
            return;
        }

        default:
        {
584
            FatalErrorInFunction
585
586
587
                << "problem." << abort(FatalError);
        }
    }
588
}
589
590


591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
void Foam::mappedPatchBase::findSamples
(
    const sampleMode mode,
    const label myWorld,
    const pointField& samples,
    const labelList& wantedWorlds,
    const labelList& origProcs,

    labelList& sampleProcs,
    labelList& sampleIndices,
    pointField& sampleLocations
) const
{
    // Find the processor/cell containing the samples. Does not account
    // for samples being found in two processors.
606

607

608
609
610
611
612
613
614
615
616
617
618
619
620
    const label myRank = Pstream::myProcNo(comm_);
    const label nProcs = Pstream::nProcs(comm_);

    wordList samplePatches(nProcs);
    {
        const label oldComm(Pstream::warnComm);
        Pstream::warnComm = comm_;
        samplePatches[myRank] = samplePatch_;
        Pstream::gatherList(samplePatches, Pstream::msgType(), comm_);
        Pstream::scatterList(samplePatches, Pstream::msgType(), comm_);
        Pstream::warnComm = oldComm;
    }
    wordList sampleRegions(nProcs);
621
    {
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
        const label oldComm(Pstream::warnComm);
        Pstream::warnComm = comm_;
        sampleRegions[myRank] = sampleRegion_;
        Pstream::gatherList(sampleRegions, Pstream::msgType(), comm_);
        Pstream::scatterList(sampleRegions, Pstream::msgType(), comm_);
        Pstream::warnComm = oldComm;
    }

    // Find all the info for nearest
    List<nearInfoWorld> nearest(samples.size());
    forAll(nearest, samplei)
    {
        nearest[samplei].first() = nearInfo
        (
            pointIndexHit(),
            Tuple2<scalar, label>(Foam::sqr(GREAT), -1)
        );
        nearest[samplei].second() = wantedWorlds[samplei];
    }

642

643
644
645
646
    // Extract samples to search for locally
    {
        DynamicList<label> localMap(samples.size());
        forAll(wantedWorlds, samplei)
647
        {
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
            if (wantedWorlds[samplei] == myWorld)
            {
                localMap.append(samplei);
            }
        }

        if (localMap.size())
        {
            pointField localSamples(samples, localMap);
            labelList localOrigProcs(origProcs, localMap);

            //Assume single patch to sample for now
            const word localOrigPatch(samplePatches[localOrigProcs[0]]);
            const word localOrigRegion(sampleRegions[localOrigProcs[0]]);
            List<nearInfoWorld> localNearest(localSamples.size());

            if (debug)
            {
                Pout<< "Searching locally for " << localSamples.size()
                    << " samples on region:" << localOrigRegion
                    << " on patch:" << localOrigPatch << endl;
            }
            findLocalSamples
            (
                mode,
                myWorld,
                localOrigRegion,
                localOrigPatch,
                localSamples,
                localNearest
            );
            UIndirectList<nearInfoWorld>(nearest, localMap) = localNearest;
680
681
682
        }
    }

683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
    const label oldComm(Pstream::warnComm);
    Pstream::warnComm = comm_;

    // Find nearest. Combine on master.
    Pstream::listCombineGather
    (
        nearest,
        nearestWorldEqOp(),
        Pstream::msgType(),
        comm_
    );
    Pstream::listCombineScatter(nearest, Pstream::msgType(), comm_);

    //if (debug)
    //{
    //    Pout<< "** AFter combining:" << endl;
    //    forAll(nearest, samplei)
    //    {
    //        Pout<< "  sample:" << samples[samplei]
    //            << " origating from proc:" << origProcs[samplei]
    //            << " to be found on world:" << wantedWorlds[samplei] << nl
    //            << "    found on world:" << nearest[samplei].second() << nl
    //            << "    found on proc:"
    //            << nearest[samplei].first().second().second() << nl
    //            << "    found on patchfacei:"
    //            << nearest[samplei].first().first().index() << nl
    //            << "    found at location:"
    //            << nearest[samplei].first().first().rawPoint() << nl;
    //    }
    //    Pout<< endl;
    //}

715
716
717
718
719
720
721
    // Convert back into proc+local index
    sampleProcs.setSize(samples.size());
    sampleIndices.setSize(samples.size());
    sampleLocations.setSize(samples.size());

    forAll(nearest, sampleI)
    {
722
723
724
        const nearInfo& ni = nearest[sampleI].first();

        if (!ni.first().hit())
725
726
727
728
729
730
731
        {
            sampleProcs[sampleI] = -1;
            sampleIndices[sampleI] = -1;
            sampleLocations[sampleI] = vector::max;
        }
        else
        {
732
733
734
            sampleProcs[sampleI] = ni.second().second();
            sampleIndices[sampleI] = ni.first().index();
            sampleLocations[sampleI] = ni.first().hitPoint();
735
        }
736
    }
737
738

    Pstream::warnComm = oldComm;
739
740
741
742
743
}


void Foam::mappedPatchBase::calcMapping() const
{
744
    static bool hasWarned = false;
745
    if (mapPtr_)
746
    {
747
        FatalErrorInFunction
748
749
750
            << "Mapping already calculated" << exit(FatalError);
    }

751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
    //// Make sure if running in database that there is a syncObjects FO
    //if (sampleDatabase() && !sameWorld())
    //{
    //    const word syncName("syncObjects");
    //    const polyMesh& mesh = patch_.boundaryMesh().mesh();
    //    const Time& runTime = mesh.time();
    //    const functionObjectList& fos = runTime.functionObjects();
    //
    //    forAll(fos, i)
    //    {
    //        Pout<< "** FO:" << fos[i].name() << " tpye:" << fos[i].type()
    //            << endl;
    //    }
    //
    //    if (fos.findObjectID(syncName) == -1)
    //    {
    //        //FatalErrorInFunction
    //        WarningInFunction
    //            << "Did not detect functionObject " << syncName
    //            << ". This is used to synchronise data inbetween worlds"
    //            << endl;
    //    }
    //}


776
777
778
779
780
    // Get points on face (since cannot use face-centres - might be off
    // face-diagonal decomposed tets.
    tmp<pointField> patchPoints(facePoints(patch_));

    // Get offsetted points
781
    const pointField offsettedPoints(samplePoints(patchPoints()));
782

783
784
    // Do a sanity check - am I sampling my own patch?
    // This only makes sense for a non-zero offset.
785
786
787
    bool sampleMyself =
    (
        mode_ == NEARESTPATCHFACE
788
     && sampleWorld() == UPstream::myWorld()
789
790
     && sampleRegion() == patch_.boundaryMesh().mesh().name()
     && samplePatch() == patch_.name()
791
792
    );

793
    if (sampleMyself)
794
    {
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
        // Check offset
        vectorField d(offsettedPoints-patchPoints());
        bool coincident = (gAverage(mag(d)) <= ROOTVSMALL);

        if (sampleMyself && coincident)
        {
            WarningInFunction
                << "Invalid offset " << d << endl
                << "Offset is the vector added to the patch face centres to"
                << " find the patch face supplying the data." << endl
                << "Setting it to " << d
                << " on the same patch, on the same region"
                << " will find the faces themselves which does not make sense"
                << " for anything but testing." << endl
                << "patch_:" << patch_.name() << endl
                << "sampleRegion_:" << sampleRegion() << endl
                << "mode_:" << sampleModeNames_[mode_] << endl
                << "samplePatch_:" << samplePatch() << endl
                << "offsetMode_:" << offsetModeNames_[offsetMode_] << endl;
        }
815
816
    }

817
818
819
820
821
    // Get local world and world-to-sample in index form
    const label myWorld = Pstream::myWorldID();
    const label mySampleWorld = Pstream::allWorlds().find(sampleWorld_);


822
    // Get global list of all samples and the processor and face they come from.
823
824
825
826
827
    pointField samples;         // coordinates
    labelList patchFaceWorlds;  // world to sample
    labelList patchFaceProcs;   // originating processor
    labelList patchFaces;       // originating face on processor patch
    pointField patchFc;         // originating face centre
828
829
    collectSamples
    (
830
        mySampleWorld,          // world I want to sample
831
        patchPoints,
832

833
        samples,
834
        patchFaceWorlds,
835
836
837
838
839
        patchFaceProcs,
        patchFaces,
        patchFc
    );

840
841
842
843
844
845
846
847
848
849
850
    //if (debug)
    //{
    //    forAll(samples, samplei)
    //    {
    //        Pout<< "    sample:" << samples[samplei]
    //            << " origating from proc:" << patchFaceProcs[samplei]
    //            << "  face:" << patchFaces[samplei]
    //            << " to be found on world:" << patchFaceWorlds[samplei] << nl;
    //    }
    //}

851
852
853
854
    // Find processor and cell/face samples are in and actual location.
    labelList sampleProcs;
    labelList sampleIndices;
    pointField sampleLocations;
855
856
857
858
859
860
861
862
863
864
865
866
    findSamples
    (
        mode_,
        myWorld,        // my world (in index form)
        samples,
        patchFaceWorlds,
        patchFaceProcs,

        sampleProcs,
        sampleIndices,
        sampleLocations
    );
867

868
869
870
871
872
873
874
875
876
877
878
879
    // Check for samples that were not found. This will only happen for
    // NEARESTCELL since finds cell containing a location
    if (mode_ == NEARESTCELL)
    {
        label nNotFound = 0;
        forAll(sampleProcs, sampleI)
        {
            if (sampleProcs[sampleI] == -1)
            {
                nNotFound++;
            }
        }
880
        reduce(nNotFound, sumOp<label>(), Pstream::msgType(), comm_);
881
882
883
884
885

        if (nNotFound > 0)
        {
            if (!hasWarned)
            {
886
887
                WarningInFunction
                    << "Did not find " << nNotFound
888
889
                    << " out of " << sampleProcs.size() << " total samples."
                    << " Sampling these on owner cell centre instead." << endl
890
                    << "On patch " << patch_.name()
891
                    << " on region " << sampleRegion()
892
                    << " in mode " << sampleModeNames_[mode_] << endl
893
894
                    << "with offset mode " << offsetModeNames_[offsetMode_]
                    << ". Suppressing further warnings from " << type() << endl;
895
896
897
898

                hasWarned = true;
            }

899
900
            // Collect the samples that cannot be found
            DynamicList<label> subMap;
901
902
903
904
            forAll(sampleProcs, sampleI)
            {
                if (sampleProcs[sampleI] == -1)
                {
905
                    subMap.append(sampleI);
906
907
908
                }
            }

909
910
911
912
913
914
915
            // And re-search for pure nearest (should not fail)
            labelList subSampleProcs;
            labelList subSampleIndices;
            pointField subSampleLocations;
            findSamples
            (
                NEARESTONLYCELL,
916
917
918
919
920
921
                myWorld,        // my world (in index form)

                pointField(samples, subMap),
                UIndirectList<label>(patchFaceWorlds, subMap)(),
                UIndirectList<label>(patchFaceProcs, subMap)(),

922
923
924
925
926
927
                subSampleProcs,
                subSampleIndices,
                subSampleLocations
            );

            // Insert
928
929
            labelUIndList(sampleProcs, subMap) = subSampleProcs;
            labelUIndList(sampleIndices, subMap) = subSampleIndices;
930
            UIndirectList<point>(sampleLocations, subMap) = subSampleLocations;
931
932
933
        }
    }

934
935
936
937
938
939
    // Now we have all the data we need:
    // - where sample originates from (so destination when mapping):
    //   patchFaces, patchFaceProcs.
    // - cell/face sample is in (so source when mapping)
    //   sampleIndices, sampleProcs.

940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
    if (Pstream::master(comm_))
    {
        forAll(samples, i)
        {
            if (sampleProcs[i] == -1)
            {
                FatalErrorInFunction << "did not find sample "
                    << patchFc[i] << " on patch " << patch_.name()
                    << " on region "
                    << patch_.boundaryMesh().mesh().name()
                    << " on processor " << patchFaceProcs[i]
                    << exit(FatalError);
            }
        }
    }
955
956


957
    if (debug && Pstream::master(comm_))
958
    {
959
960
961
962
963
964
965
966
967
968
969
970
971
972
        //forAll(samples, i)
        //{
        //    Pout<< i << " need data in region "
        //        << patch_.boundaryMesh().mesh().name()
        //        << " for proc:" << patchFaceProcs[i]
        //        << " face:" << patchFaces[i]
        //        << " at:" << patchFc[i] << endl
        //        << "Found data in region " << sampleRegion()
        //        << " at proc:" << sampleProcs[i]
        //        << " face:" << sampleIndices[i]
        //        << " at:" << sampleLocations[i]
        //        << nl << endl;
        //}

973
974
975
976
977
978
979
        OFstream str
        (
            patch_.boundaryMesh().mesh().time().path()
          / patch_.name()
          + "_mapped.obj"
        );
        Pout<< "Dumping mapping as lines from patch faceCentres to"
980
981
            << " sampled cell/faceCentres/points to file " << str.name()
            << endl;
982
983
984
985
986
987
988
989
990
991
992
993
994
995

        label vertI = 0;

        forAll(patchFc, i)
        {
            meshTools::writeOBJ(str, patchFc[i]);
            vertI++;
            meshTools::writeOBJ(str, sampleLocations[i]);
            vertI++;
            str << "l " << vertI-1 << ' ' << vertI << nl;
        }
    }

    // Determine schedule.
996
    mapPtr_.reset(new mapDistribute(sampleProcs, patchFaceProcs, comm_));
997
998
999
1000
1001
1002
1003

    // Rework the schedule from indices into samples to cell data to send,
    // face data to receive.

    labelListList& subMap = mapPtr_().subMap();
    labelListList& constructMap = mapPtr_().constructMap();

1004
    forAll(subMap, proci)
1005
    {
1006
1007
        subMap[proci] = labelUIndList(sampleIndices, subMap[proci]);
        constructMap[proci] = labelUIndList(patchFaces, constructMap[proci]);
1008

1009
1010
1011
1012
1013
1014
1015
1016
        if (debug)
        {
            Pout<< "To proc:" << proci << " sending values of cells/faces:"
                << subMap[proci] << endl;
            Pout<< "From proc:" << proci
                << " receiving values of patch faces:"
                << constructMap[proci] << endl;
        }
1017
1018
    }

1019

1020
1021
1022
1023
1024
1025
    // Redo constructSize
    mapPtr_().constructSize() = patch_.size();

    if (debug)
    {
        // Check that all elements get a value.
1026
        bitSet used(patch_.size());
1027
        forAll(constructMap, proci)
1028
        {
1029
            const labelList& map = constructMap[proci];
1030
1031
1032

            forAll(map, i)
            {
1033
                label facei = map[i];
1034

1035
                if (used.test(facei))
1036
                {
1037
                    FatalErrorInFunction
1038
                        << "On patch " << patch_.name()
1039
                        << " patchface " << facei
1040
1041
1042
                        << " is assigned to more than once."
                        << abort(FatalError);
                }
1043
1044
1045
1046
                else
                {
                    used.set(facei);
                }
1047
1048
            }
        }
1049
        forAll(used, facei)
1050
        {
1051
            if (!used.test(facei))
1052
            {
1053
                FatalErrorInFunction
1054
                    << "On patch " << patch_.name()
1055
                    << " patchface " << facei
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
                    << " is never assigned to."
                    << abort(FatalError);
            }
        }
    }
}


const Foam::autoPtr<Foam::searchableSurface>& Foam::mappedPatchBase::surfPtr()
const
{
1067
    const word surfType(surfDict_.getOrDefault<word>("type", "none"));
1068

1069
    if (!surfPtr_ && surfType != "none")
1070
    {
1071
        word surfName(surfDict_.getOrDefault("name", patch_.name()));
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097

        const polyMesh& mesh = patch_.boundaryMesh().mesh();

        surfPtr_ =
            searchableSurface::New
            (
                surfType,
                IOobject
                (
                    surfName,
                    mesh.time().constant(),
                    "triSurface",
                    mesh,
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE
                ),
                surfDict_
            );
    }

    return surfPtr_;
}


void Foam::mappedPatchBase::calcAMI() const
{
Andrew Heather's avatar
Andrew Heather committed
1098
    if (AMIPtr_->upToDate())
1099
    {
1100
        DebugInFunction
Andrew Heather's avatar
Andrew Heather committed
1101
1102
            << "AMI already up-to-date"
            << endl;
1103
1104

        return;
1105
1106
    }

1107
1108
1109
1110
    // Check if running locally
    if (sampleWorld_.empty())
    {
        const polyPatch& nbr = samplePolyPatch();
1111

1112
1113
        // Transform neighbour patch to local system
        pointField nbrPoints(samplePoints(nbr.localPoints()));
1114

1115
        primitivePatch nbrPatch0
1116
        (
1117
1118
1119
1120
1121
1122
1123
            SubList<face>
            (
                nbr.localFaces(),
                nbr.size()
            ),
            nbrPoints
        );
1124
1125


1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
        if (debug)
        {
            OFstream os(patch_.name() + "_neighbourPatch-org.obj");
            meshTools::writeOBJ(os, samplePolyPatch().localFaces(), nbrPoints);

            OFstream osN(patch_.name() + "_neighbourPatch-trans.obj");
            meshTools::writeOBJ(osN, nbrPatch0, nbrPoints);

            OFstream osO(patch_.name() + "_ownerPatch.obj");
            meshTools::writeOBJ(osO, patch_.localFaces(), patch_.localPoints());
        }

        // Construct/apply AMI interpolation to determine addressing and weights
        AMIPtr_->calculate(patch_, nbrPatch0, surfPtr());
    }
    else
1142
    {
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
        faceList nbrFaces;
        pointField nbrPoints;
        if (sampleWorld() == UPstream::myWorld())
        {
            const polyPatch& nbr = samplePolyPatch();
            nbrFaces = nbr.localFaces();
            nbrPoints = samplePoints(nbr.localPoints());
        }
        else
        {
            // Leave empty
        }
1155

1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
        primitivePatch nbrPatch0
        (
            SubList<face>
            (
                nbrFaces,
                nbrFaces.size()
            ),
            nbrPoints
        );

        // Change to use ALL processors communicator
        const label oldWorldComm = Pstream::worldComm;
        Pstream::worldComm = comm_;

        const label oldComm(Pstream::warnComm);
        Pstream::warnComm = UPstream::worldComm;
1172

1173
1174
1175
1176
1177
        // Construct/apply AMI interpolation to determine addressing and weights
        AMIPtr_->calculate(patch_, nbrPatch0, surfPtr());

        Pstream::warnComm = oldComm;
        Pstream::worldComm = oldWorldComm;
1178
    }
1179
}
1180

1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197

const Foam::objectRegistry& Foam::mappedPatchBase::subRegistry
(
    const objectRegistry& obr,
    const wordList& names,
    const label index
)
{
    const objectRegistry& sub = obr.subRegistry(names[index], true);
    if (index == names.size()-1)
    {
        return sub;
    }
    else
    {
        return subRegistry(sub, names, index+1);
    }
1198
1199
1200
1201
1202
}


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

Andrew Heather's avatar
Andrew Heather committed
1203
Foam::mappedPatchBase::mappedPatchBase(const polyPatch& pp)
1204
1205
:
    patch_(pp),
1206
    sampleWorld_(word::null),
1207
1208
    sampleRegion_(patch_.boundaryMesh().mesh().name()),
    mode_(NEARESTPATCHFACE),
1209
    samplePatch_(word::null),
1210
    coupleGroup_(),
1211
    sampleDatabasePtr_(),
1212
    offsetMode_(UNIFORM),
Henry Weller's avatar
Henry Weller committed
1213
    offset_(Zero),
1214
1215
    offsets_(pp.size(), offset_),
    distance_(0),
1216
1217
1218
1219
1220
1221
    comm_(communicator(sampleWorld_)),
    sameRegion_
    (
        sampleWorld_.empty()
     && sampleRegion_ == patch_.boundaryMesh().mesh().name()
    ),
1222
    mapPtr_(nullptr),
1223
    AMIReverse_(false),
Andrew Heather's avatar
Andrew Heather committed
1224
    AMIPtr_(new faceAreaWeightAMI(true, AMIReverse_)),
1225
    surfPtr_(nullptr),
1226
    surfDict_(fileName("surface"))
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
    const word& sampleRegion,
    const sampleMode mode,
    const word& samplePatch,
    const vectorField& offsets
)
:
    patch_(pp),
1240
    sampleWorld_(word::null),
1241
1242
1243
    sampleRegion_(sampleRegion),
    mode_(mode),
    samplePatch_(samplePatch),
1244
    coupleGroup_(),
1245
    sampleDatabasePtr_(),
1246
    offsetMode_(NONUNIFORM),
Henry Weller's avatar
Henry Weller committed
1247
    offset_(Zero),
1248
1249
    offsets_(offsets),
    distance_(0),
1250
1251
1252
1253
1254
1255
    comm_(communicator(sampleWorld_)),
    sameRegion_
    (
        sampleWorld_.empty()
     && sampleRegion_ == patch_.boundaryMesh().mesh().name()
    ),
1256
    mapPtr_(nullptr),
1257
    AMIReverse_(false),
Andrew Heather's avatar
Andrew Heather committed
1258
    AMIPtr_(new faceAreaWeightAMI(true, AMIReverse_)),
1259
    surfPtr_(nullptr),
1260
    surfDict_(fileName("surface"))
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
    const word& sampleRegion,
    const sampleMode mode,
    const word& samplePatch,
    const vector& offset
)
:
    patch_(pp),
1274
    sampleWorld_(word::null),
1275
1276
1277
    sampleRegion_(sampleRegion),
    mode_(mode),
    samplePatch_(samplePatch),
1278
    coupleGroup_(),
1279
    sampleDatabasePtr_(),
1280
1281
1282
1283
    offsetMode_(UNIFORM),
    offset_(offset),
    offsets_(0),
    distance_(0),
1284
1285
1286
1287
1288
1289
    comm_(communicator(sampleWorld_)),
    sameRegion_
    (
        sampleWorld_.empty()
     && sampleRegion_ == patch_.boundaryMesh().mesh().name()
    ),
1290
    mapPtr_(nullptr),
1291
    AMIReverse_(false),
Andrew Heather's avatar
Andrew Heather committed
1292
    AMIPtr_(new faceAreaWeightAMI(true, AMIReverse_)),
1293
    surfPtr_(nullptr),
1294
    surfDict_(fileName("surface"))
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
    const word& sampleRegion,
    const sampleMode mode,
    const word& samplePatch,
    const scalar distance
)
:
    patch_(pp),
1308
    sampleWorld_(word::null),
1309
1310
1311
    sampleRegion_(sampleRegion),
    mode_(mode),
    samplePatch_(samplePatch),
1312
    coupleGroup_(),
1313
    sampleDatabasePtr_(),
1314
    offsetMode_(NORMAL),
Henry Weller's avatar
Henry Weller committed
1315
    offset_(Zero),
1316
1317
    offsets_(0),
    distance_(distance),
1318
1319
1320
1321
1322
1323
    comm_(communicator(sampleWorld_)),
    sameRegion_
    (
        sampleWorld_.empty()
     && sampleRegion_ == patch_.boundaryMesh().mesh().name()
    ),
1324
    mapPtr_(nullptr),
1325
    AMIReverse_(false),
Andrew Heather's avatar
Andrew Heather committed
1326
    AMIPtr_(new faceAreaWeightAMI(true, AMIReverse_)),
1327
    surfPtr_(nullptr),
1328
    surfDict_(fileName("surface"))
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
    const dictionary& dict
)
:
    patch_(pp),
1339
1340
    sampleWorld_(dict.getOrDefault<word>("sampleWorld", word::null)),
    sampleRegion_(dict.getOrDefault<word>("sampleRegion", word::null)),
Mark OLESEN's avatar
Mark OLESEN committed
1341
    mode_(sampleModeNames_.get("sampleMode", dict)),
1342
    samplePatch_(dict.getOrDefault<word>("samplePatch", word::null)),
1343
    coupleGroup_(dict),
1344
    sampleDatabasePtr_(readDatabase(dict)),
1345
    offsetMode_(UNIFORM),