mappedPatchBase.C 39.5 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
6
     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
-------------------------------------------------------------------------------
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"
29
#include "meshSearchMeshObject.H"
30
31
32
33
#include "meshTools.H"
#include "OFstream.H"
#include "Random.H"
#include "treeDataFace.H"
34
#include "treeDataPoint.H"
35
36
37
38
39
40
#include "indexedOctree.H"
#include "polyMesh.H"
#include "polyPatch.H"
#include "Time.H"
#include "mapDistribute.H"
#include "SubField.H"
41
#include "triPointRef.H"
42
#include "syncTools.H"
43
#include "treeDataCell.H"
44
#include "DynamicField.H"
45
46
47
48
49
50
51
52
53

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

namespace Foam
{
    defineTypeNameAndDebug(mappedPatchBase, 0);
}


54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
const Foam::Enum
<
    Foam::mappedPatchBase::sampleMode
>
Foam::mappedPatchBase::sampleModeNames_
{
    { sampleMode::NEARESTCELL, "nearestCell" },
    { sampleMode::NEARESTPATCHFACE, "nearestPatchFace" },
    { sampleMode::NEARESTPATCHFACEAMI, "nearestPatchFaceAMI" },
    { sampleMode::NEARESTPATCHPOINT, "nearestPatchPoint" },
    { sampleMode::NEARESTFACE, "nearestFace" },
    { sampleMode::NEARESTONLYCELL, "nearestOnlyCell" },
};


const Foam::Enum
<
    Foam::mappedPatchBase::offsetMode
>
Foam::mappedPatchBase::offsetModeNames_
{
    { offsetMode::UNIFORM, "uniform" },
    { offsetMode::NONUNIFORM, "nonuniform" },
    { offsetMode::NORMAL, "normal" },
};
79
80
81
82


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

83
84
85
86
87
88
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::facePoints
(
    const polyPatch& pp
) const
{
    const polyMesh& mesh = pp.boundaryMesh().mesh();
89
90
91

    // Force construction of min-tet decomp
    (void)mesh.tetBasePtIs();
92
93

    // Initialise to face-centre
94
    tmp<pointField> tfacePoints(new pointField(patch_.size()));
95
    pointField& facePoints = tfacePoints.ref();
96

97
    forAll(pp, facei)
98
    {
99
        facePoints[facei] = facePoint
100
101
        (
            mesh,
102
            pp.start()+facei,
103
            polyMesh::FACE_DIAG_TRIS
104
        ).rawPoint();
105
106
107
108
109
110
    }

    return tfacePoints;
}


111
112
void Foam::mappedPatchBase::collectSamples
(
113
    const pointField& facePoints,
114
115
116
117
118
119
120
    pointField& samples,
    labelList& patchFaceProcs,
    labelList& patchFaces,
    pointField& patchFc
) const
{
    // Collect all sample points and the faces they come from.
121
122
123
124
125
126
127
128
129
130
131
132
    {
        List<pointField> globalFc(Pstream::nProcs());
        globalFc[Pstream::myProcNo()] = facePoints;
        Pstream::gatherList(globalFc);
        Pstream::scatterList(globalFc);
        // Rework into straight list
        patchFc = ListListOps::combine<pointField>
        (
            globalFc,
            accessOp<pointField>()
        );
    }
133

134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
    {
        List<pointField> globalSamples(Pstream::nProcs());
        globalSamples[Pstream::myProcNo()] = samplePoints(facePoints);
        Pstream::gatherList(globalSamples);
        Pstream::scatterList(globalSamples);
        // Rework into straight list
        samples = ListListOps::combine<pointField>
        (
            globalSamples,
            accessOp<pointField>()
        );
    }

    {
        labelListList globalFaces(Pstream::nProcs());
        globalFaces[Pstream::myProcNo()] = identity(patch_.size());
        // Distribute to all processors
        Pstream::gatherList(globalFaces);
        Pstream::scatterList(globalFaces);

        patchFaces = ListListOps::combine<labelList>
155
156
157
        (
            globalFaces,
            accessOp<labelList>()
158
159
160
        );
    }

161
    {
162
163
164
165
166
167
168
169
        labelList nPerProc(Pstream::nProcs());
        nPerProc[Pstream::myProcNo()] = patch_.size();
        Pstream::gatherList(nPerProc);
        Pstream::scatterList(nPerProc);

        patchFaceProcs.setSize(patchFaces.size());

        label sampleI = 0;
170
        forAll(nPerProc, proci)
171
        {
172
            for (label i = 0; i < nPerProc[proci]; i++)
173
            {
174
                patchFaceProcs[sampleI++] = proci;
175
            }
176
177
178
179
180
181
182
183
184
        }
    }
}


// Find the processor/cell containing the samples. Does not account
// for samples being found in two processors.
void Foam::mappedPatchBase::findSamples
(
185
    const sampleMode mode,
186
187
188
189
190
191
192
193
194
195
196
197
    const pointField& samples,
    labelList& sampleProcs,
    labelList& sampleIndices,
    pointField& sampleLocations
) const
{
    // Lookup the correct region
    const polyMesh& mesh = sampleMesh();

    // All the info for nearest. Construct to miss
    List<nearInfo> nearest(samples.size());

198
    switch (mode)
199
200
201
    {
        case NEARESTCELL:
        {
202
            if (samplePatch_.size() && samplePatch_ != "none")
203
            {
204
205
                FatalErrorInFunction
                    << "No need to supply a patch name when in "
206
                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
207
208
            }

209
            //- Note: face-diagonal decomposition
210
            const indexedOctree<Foam::treeDataCell>& tree = mesh.cellTree();
211
212
213
214
215

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

216
                label celli = tree.findInside(sample);
217

218
                if (celli == -1)
219
220
221
222
223
224
                {
                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
                    nearest[sampleI].second().second() = Pstream::myProcNo();
                }
                else
                {
225
                    const point& cc = mesh.cellCentres()[celli];
226
227
228
229
230

                    nearest[sampleI].first() = pointIndexHit
                    (
                        true,
                        cc,
231
                        celli
232
233
234
235
236
237
238
239
                    );
                    nearest[sampleI].second().first() = magSqr(cc-sample);
                    nearest[sampleI].second().second() = Pstream::myProcNo();
                }
            }
            break;
        }

240
241
242
243
        case NEARESTONLYCELL:
        {
            if (samplePatch_.size() && samplePatch_ != "none")
            {
244
245
                FatalErrorInFunction
                    << "No need to supply a patch name when in "
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
            }

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

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

                nearest[sampleI].first() = tree.findNearest(sample, sqr(GREAT));
                nearest[sampleI].second().first() = magSqr
                (
                    nearest[sampleI].first().hitPoint()
                   -sample
                );
                nearest[sampleI].second().second() = Pstream::myProcNo();
            }
            break;
        }

267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
        case NEARESTPATCHFACE:
        {
            Random rndGen(123456);

            const polyPatch& pp = samplePolyPatch();

            if (pp.empty())
            {
                forAll(samples, sampleI)
                {
                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
                    nearest[sampleI].second().second() = Pstream::myProcNo();
                }
            }
            else
            {
                // patch faces
284
                const labelList patchFaces(identity(pp.size(), pp.start()));
285
286
287
288
289
290

                treeBoundBox patchBb
                (
                    treeBoundBox(pp.points(), pp.meshPoints()).extend
                    (
                        rndGen,
291
                        1e-4
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
                    )
                );
                patchBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
                patchBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);

                indexedOctree<treeDataFace> boundaryTree
                (
                    treeDataFace    // all information needed to search faces
                    (
                        false,      // do not cache bb
                        mesh,
                        patchFaces  // boundary faces only
                    ),
                    patchBb,        // overall search domain
                    8,              // maxLevel
                    10,             // leafsize
                    3.0             // duplicity
                );

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

                    pointIndexHit& nearInfo = nearest[sampleI].first();
                    nearInfo = boundaryTree.findNearest
                    (
                        sample,
                        magSqr(patchBb.span())
                    );

                    if (!nearInfo.hit())
                    {
                        nearest[sampleI].second().first() = Foam::sqr(GREAT);
                        nearest[sampleI].second().second() =
                            Pstream::myProcNo();
                    }
                    else
                    {
                        point fc(pp[nearInfo.index()].centre(pp.points()));
                        nearInfo.setPoint(fc);
                        nearest[sampleI].second().first() = magSqr(fc-sample);
                        nearest[sampleI].second().second() =
                            Pstream::myProcNo();
                    }
                }
            }
            break;
        }

341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
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
386
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
        case NEARESTPATCHPOINT:
        {
            Random rndGen(123456);

            const polyPatch& pp = samplePolyPatch();

            if (pp.empty())
            {
                forAll(samples, sampleI)
                {
                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
                    nearest[sampleI].second().second() = Pstream::myProcNo();
                }
            }
            else
            {
                // patch (local) points
                treeBoundBox patchBb
                (
                    treeBoundBox(pp.points(), pp.meshPoints()).extend
                    (
                        rndGen,
                        1e-4
                    )
                );
                patchBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
                patchBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);

                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];

                    pointIndexHit& nearInfo = nearest[sampleI].first();
                    nearInfo = boundaryTree.findNearest
                    (
                        sample,
                        magSqr(patchBb.span())
                    );

                    if (!nearInfo.hit())
                    {
                        nearest[sampleI].second().first() = Foam::sqr(GREAT);
                        nearest[sampleI].second().second() =
                            Pstream::myProcNo();
                    }
                    else
                    {
                        const point& pt = nearInfo.hitPoint();

                        nearest[sampleI].second().first() = magSqr(pt-sample);
                        nearest[sampleI].second().second() =
                            Pstream::myProcNo();
                    }
                }
            }
            break;
        }

412
413
        case NEARESTFACE:
        {
414
            if (samplePatch().size() && samplePatch() != "none")
415
            {
416
417
                FatalErrorInFunction
                    << "No need to supply a patch name when in "
418
                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
419
420
            }

421
422
423
            //- Note: face-diagonal decomposition
            const meshSearchMeshObject& meshSearchEngine =
                meshSearchMeshObject::New(mesh);
424
425
426
427
428

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

429
                label facei = meshSearchEngine.findNearestFace(sample);
430

431
                if (facei == -1)
432
433
434
435
436
437
                {
                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
                    nearest[sampleI].second().second() = Pstream::myProcNo();
                }
                else
                {
438
                    const point& fc = mesh.faceCentres()[facei];
439
440
441
442
443

                    nearest[sampleI].first() = pointIndexHit
                    (
                        true,
                        fc,
444
                        facei
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
                    );
                    nearest[sampleI].second().first() = magSqr(fc-sample);
                    nearest[sampleI].second().second() = Pstream::myProcNo();
                }
            }
            break;
        }

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

        default:
        {
461
            FatalErrorInFunction
462
463
464
465
466
                << "problem." << abort(FatalError);
        }
    }


467
    // Find nearest. Combine on master.
468
469
470
    Pstream::listCombineGather(nearest, nearestEqOp());
    Pstream::listCombineScatter(nearest);

471

472
473
    if (debug)
    {
474
475
476
        InfoInFunction
            << "mesh " << sampleRegion() << " : " << endl;

477
478
        forAll(nearest, sampleI)
        {
479
            label proci = nearest[sampleI].second().second();
480
481
482
            label localI = nearest[sampleI].first().index();

            Info<< "    " << sampleI << " coord:"<< samples[sampleI]
483
                << " found on processor:" << proci
484
485
486
                << " in local cell/face/point:" << localI
                << " with location:" << nearest[sampleI].first().rawPoint()
                << endl;
487
488
489
490
491
492
493
494
495
496
        }
    }

    // Convert back into proc+local index
    sampleProcs.setSize(samples.size());
    sampleIndices.setSize(samples.size());
    sampleLocations.setSize(samples.size());

    forAll(nearest, sampleI)
    {
497
498
499
500
501
502
503
504
505
506
507
508
        if (!nearest[sampleI].first().hit())
        {
            sampleProcs[sampleI] = -1;
            sampleIndices[sampleI] = -1;
            sampleLocations[sampleI] = vector::max;
        }
        else
        {
            sampleProcs[sampleI] = nearest[sampleI].second().second();
            sampleIndices[sampleI] = nearest[sampleI].first().index();
            sampleLocations[sampleI] = nearest[sampleI].first().hitPoint();
        }
509
510
511
512
513
514
    }
}


void Foam::mappedPatchBase::calcMapping() const
{
515
    static bool hasWarned = false;
516
517
    if (mapPtr_.valid())
    {
518
        FatalErrorInFunction
519
520
521
            << "Mapping already calculated" << exit(FatalError);
    }

522
523
524
525
526
    // Get points on face (since cannot use face-centres - might be off
    // face-diagonal decomposed tets.
    tmp<pointField> patchPoints(facePoints(patch_));

    // Get offsetted points
527
    const pointField offsettedPoints(samplePoints(patchPoints()));
528

529
530
    // Do a sanity check - am I sampling my own patch?
    // This only makes sense for a non-zero offset.
531
532
533
    bool sampleMyself =
    (
        mode_ == NEARESTPATCHFACE
534
535
     && sampleRegion() == patch_.boundaryMesh().mesh().name()
     && samplePatch() == patch_.name()
536
537
538
    );

    // Check offset
539
540
541
542
    vectorField d(offsettedPoints-patchPoints());
    bool coincident = (gAverage(mag(d)) <= ROOTVSMALL);

    if (sampleMyself && coincident)
543
    {
544
545
        WarningInFunction
            << "Invalid offset " << d << endl
546
547
548
549
550
551
552
            << "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
553
            << "sampleRegion_:" << sampleRegion() << endl
554
            << "mode_:" << sampleModeNames_[mode_] << endl
555
            << "samplePatch_:" << samplePatch() << endl
556
557
558
559
560
561
562
563
            << "offsetMode_:" << offsetModeNames_[offsetMode_] << endl;
    }

    // Get global list of all samples and the processor and face they come from.
    pointField samples;
    labelList patchFaceProcs;
    labelList patchFaces;
    pointField patchFc;
564
565
566
567
568
569
570
571
572
    collectSamples
    (
        patchPoints,
        samples,
        patchFaceProcs,
        patchFaces,
        patchFc
    );

573
574
575
576
    // Find processor and cell/face samples are in and actual location.
    labelList sampleProcs;
    labelList sampleIndices;
    pointField sampleLocations;
577
    findSamples(mode_, samples, sampleProcs, sampleIndices, sampleLocations);
578

579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
    // 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++;
            }
        }
        reduce(nNotFound, sumOp<label>());

        if (nNotFound > 0)
        {
            if (!hasWarned)
            {
597
598
                WarningInFunction
                    << "Did not find " << nNotFound
599
600
                    << " out of " << sampleProcs.size() << " total samples."
                    << " Sampling these on owner cell centre instead." << endl
601
                    << "On patch " << patch_.name()
602
                    << " on region " << sampleRegion()
603
                    << " in mode " << sampleModeNames_[mode_] << endl
604
605
                    << "with offset mode " << offsetModeNames_[offsetMode_]
                    << ". Suppressing further warnings from " << type() << endl;
606
607
608
609

                hasWarned = true;
            }

610
611
612
            // Collect the samples that cannot be found
            DynamicList<label> subMap;
            DynamicField<point> subSamples;
613
614
615
616
617

            forAll(sampleProcs, sampleI)
            {
                if (sampleProcs[sampleI] == -1)
                {
618
619
                    subMap.append(sampleI);
                    subSamples.append(samples[sampleI]);
620
621
622
                }
            }

623
624
625
626
627
628
629
630
631
632
633
634
635
636
            // And re-search for pure nearest (should not fail)
            labelList subSampleProcs;
            labelList subSampleIndices;
            pointField subSampleLocations;
            findSamples
            (
                NEARESTONLYCELL,
                subSamples,
                subSampleProcs,
                subSampleIndices,
                subSampleLocations
            );

            // Insert
637
638
            labelUIndList(sampleProcs, subMap) = subSampleProcs;
            labelUIndList(sampleIndices, subMap) = subSampleIndices;
639
            UIndirectList<point>(sampleLocations, subMap) = subSampleLocations;
640
641
642
        }
    }

643
644
645
646
647
648
649
650
651
652
653
654
655
    // 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.

    //forAll(samples, i)
    //{
    //    Info<< i << " need data in region "
    //        << patch_.boundaryMesh().mesh().name()
    //        << " for proc:" << patchFaceProcs[i]
    //        << " face:" << patchFaces[i]
    //        << " at:" << patchFc[i] << endl
656
    //        << "Found data in region " << sampleRegion()
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
    //        << " at proc:" << sampleProcs[i]
    //        << " face:" << sampleIndices[i]
    //        << " at:" << sampleLocations[i]
    //        << nl << endl;
    //}



    if (debug && Pstream::master())
    {
        OFstream str
        (
            patch_.boundaryMesh().mesh().time().path()
          / patch_.name()
          + "_mapped.obj"
        );
        Pout<< "Dumping mapping as lines from patch faceCentres to"
674
675
            << " sampled cell/faceCentres/points to file " << str.name()
            << endl;
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697

        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.
    mapPtr_.reset(new mapDistribute(sampleProcs, patchFaceProcs));

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

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

698
    forAll(subMap, proci)
699
    {
700
701
        subMap[proci] = labelUIndList(sampleIndices, subMap[proci]);
        constructMap[proci] = labelUIndList(patchFaces, constructMap[proci]);
702
703
704

        //if (debug)
        //{
705
706
707
        //    Pout<< "To proc:" << proci << " sending values of cells/faces:"
        //        << subMap[proci] << endl;
        //    Pout<< "From proc:" << proci
708
        //        << " receiving values of patch faces:"
709
        //        << constructMap[proci] << endl;
710
711
712
713
714
715
716
717
718
719
        //}
    }

    // Redo constructSize
    mapPtr_().constructSize() = patch_.size();

    if (debug)
    {
        // Check that all elements get a value.
        PackedBoolList used(patch_.size());
720
        forAll(constructMap, proci)
721
        {
722
            const labelList& map = constructMap[proci];
723
724
725

            forAll(map, i)
            {
726
                label facei = map[i];
727

728
                if (used[facei] == 0)
729
                {
730
                    used[facei] = 1;
731
732
733
                }
                else
                {
734
                    FatalErrorInFunction
735
                        << "On patch " << patch_.name()
736
                        << " patchface " << facei
737
738
739
740
741
                        << " is assigned to more than once."
                        << abort(FatalError);
                }
            }
        }
742
        forAll(used, facei)
743
        {
744
            if (used[facei] == 0)
745
            {
746
                FatalErrorInFunction
747
                    << "On patch " << patch_.name()
748
                    << " patchface " << facei
749
750
751
752
753
754
755
756
757
758
759
                    << " is never assigned to."
                    << abort(FatalError);
            }
        }
    }
}


const Foam::autoPtr<Foam::searchableSurface>& Foam::mappedPatchBase::surfPtr()
const
{
760
    const word surfType(surfDict_.lookupOrDefault<word>("type", "none"));
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792

    if (!surfPtr_.valid() && surfType != "none")
    {
        word surfName(surfDict_.lookupOrDefault("name", patch_.name()));

        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
{
    if (AMIPtr_.valid())
    {
793
        FatalErrorInFunction
794
795
796
797
798
            << "AMI already calculated" << exit(FatalError);
    }

    AMIPtr_.clear();

799

800
    const polyPatch& nbr = samplePolyPatch();
801

802
803
    // Transform neighbour patch to local system
    pointField nbrPoints(samplePoints(nbr.localPoints()));
804

805
806
807
    primitivePatch nbrPatch0
    (
        SubList<face>
808
        (
809
810
811
812
813
814
815
816
817
818
819
            nbr.localFaces(),
            nbr.size()
        ),
        nbrPoints
    );


    if (debug)
    {
        OFstream os(patch_.name() + "_neighbourPatch-org.obj");
        meshTools::writeOBJ(os, samplePolyPatch().localFaces(), nbrPoints);
820
821
822
823
824
825
826

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

828
829
830
831
832
833
    // Construct/apply AMI interpolation to determine addressing and weights
    AMIPtr_.reset
    (
        new AMIPatchToPatchInterpolation
        (
            patch_,
834
            nbrPatch0,
835
            surfPtr(),
836
            faceAreaIntersect::tmMesh,
837
            true,
andy's avatar
andy committed
838
            AMIPatchToPatchInterpolation::imFaceAreaWeight,
839
            -1,
840
            AMIReverse_
841
842
843
844
845
        )
    );
}


846
847
848
849
850
851
852
853
854
855
856
// Hack to read old (List-based) format. See Field.C. The difference
// is only that in case of falling back to old format it expects a non-uniform
// list instead of a single vector.
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::readListOrField
(
    const word& keyword,
    const dictionary& dict,
    const label size
)
{
    tmp<pointField> tfld(new pointField());
857
    pointField& fld = tfld.ref();
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877

    if (size)
    {
        ITstream& is = dict.lookup(keyword);

        // Read first token
        token firstToken(is);

        if (firstToken.isWord())
        {
            if (firstToken.wordToken() == "uniform")
            {
                fld.setSize(size);
                fld = pTraits<vector>(is);
            }
            else if (firstToken.wordToken() == "nonuniform")
            {
                is >> static_cast<List<vector>&>(fld);
                if (fld.size() != size)
                {
878
                    FatalIOErrorInFunction
879
880
881
882
883
884
885
886
887
                    (
                        dict
                    )   << "size " << fld.size()
                        << " is not equal to the given value of " << size
                        << exit(FatalIOError);
                }
            }
            else
            {
888
                FatalIOErrorInFunction
889
890
891
892
893
894
895
896
897
898
899
                (
                    dict
                )   << "expected keyword 'uniform' or 'nonuniform', found "
                    << firstToken.wordToken()
                    << exit(FatalIOError);
            }
        }
        else
        {
            if (is.version() == 2.0)
            {
900
                IOWarningInFunction
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
                (
                    dict
                )   << "expected keyword 'uniform' or 'nonuniform', "
                       "assuming List format for backwards compatibility."
                       "Foam version 2.0." << endl;

                is.putBack(firstToken);
                is >> static_cast<List<vector>&>(fld);
            }
        }
    }
    return tfld;
}


916
917
918
919
920
921
922
923
924
925
// * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //

Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp
)
:
    patch_(pp),
    sampleRegion_(patch_.boundaryMesh().mesh().name()),
    mode_(NEARESTPATCHFACE),
926
927
    samplePatch_(""),
    coupleGroup_(),
928
    offsetMode_(UNIFORM),
Henry Weller's avatar
Henry Weller committed
929
    offset_(Zero),
930
931
932
    offsets_(pp.size(), offset_),
    distance_(0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
933
934
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
935
    AMIReverse_(false),
936
    surfPtr_(nullptr),
937
    surfDict_(fileName("surface"))
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
    const word& sampleRegion,
    const sampleMode mode,
    const word& samplePatch,
    const vectorField& offsets
)
:
    patch_(pp),
    sampleRegion_(sampleRegion),
    mode_(mode),
    samplePatch_(samplePatch),
954
    coupleGroup_(),
955
    offsetMode_(NONUNIFORM),
Henry Weller's avatar
Henry Weller committed
956
    offset_(Zero),
957
958
959
    offsets_(offsets),
    distance_(0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
960
961
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
962
    AMIReverse_(false),
963
    surfPtr_(nullptr),
964
    surfDict_(fileName("surface"))
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
    const word& sampleRegion,
    const sampleMode mode,
    const word& samplePatch,
    const vector& offset
)
:
    patch_(pp),
    sampleRegion_(sampleRegion),
    mode_(mode),
    samplePatch_(samplePatch),
981
    coupleGroup_(),
982
983
984
985
986
    offsetMode_(UNIFORM),
    offset_(offset),
    offsets_(0),
    distance_(0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
987
988
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
989
    AMIReverse_(false),
990
    surfPtr_(nullptr),
991
    surfDict_(fileName("surface"))
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
    const word& sampleRegion,
    const sampleMode mode,
    const word& samplePatch,
    const scalar distance
)
:
    patch_(pp),
    sampleRegion_(sampleRegion),
    mode_(mode),
    samplePatch_(samplePatch),
1008
    coupleGroup_(),
1009
    offsetMode_(NORMAL),
Henry Weller's avatar
Henry Weller committed
1010
    offset_(Zero),
1011
1012
1013
    offsets_(0),
    distance_(distance),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
1014
1015
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1016
    AMIReverse_(false),
1017
    surfPtr_(nullptr),
1018
    surfDict_(fileName("surface"))
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
    const dictionary& dict
)
:
    patch_(pp),
1029
    sampleRegion_(dict.lookupOrDefault<word>("sampleRegion", "")),
1030
    mode_(sampleModeNames_.lookup("sampleMode", dict)),
1031
1032
    samplePatch_(dict.lookupOrDefault<word>("samplePatch", "")),
    coupleGroup_(dict),
1033
    offsetMode_(UNIFORM),
Henry Weller's avatar
Henry Weller committed
1034
    offset_(Zero),
1035
1036
1037
    offsets_(0),
    distance_(0.0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
1038
1039
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1040
    AMIReverse_(dict.lookupOrDefault<bool>("flipNormals", false)),
1041
    surfPtr_(nullptr),
1042
1043
    surfDict_(dict.subOrEmptyDict("surface"))
{
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
    if (!coupleGroup_.valid())
    {
        if (sampleRegion_.empty())
        {
            // If no coupleGroup and no sampleRegion assume local region
            sampleRegion_ = patch_.boundaryMesh().mesh().name();
            sameRegion_ = true;
        }
    }

1054
1055
    if (dict.found("offsetMode"))
    {
1056
        offsetMode_ = offsetModeNames_.lookup("offsetMode", dict);
1057

Mark OLESEN's avatar
Mark OLESEN committed
1058
        switch (offsetMode_)
1059
1060
1061
1062
1063
1064
1065
1066
1067
        {
            case UNIFORM:
            {
                offset_ = point(dict.lookup("offset"));
            }
            break;

            case NONUNIFORM:
            {
1068
1069
                //offsets_ = pointField(dict.lookup("offsets"));
                offsets_ = readListOrField("offsets", dict, patch_.size());
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
            }
            break;

            case NORMAL:
            {
                distance_ = readScalar(dict.lookup("distance"));
            }
            break;
        }
    }
    else if (dict.found("offset"))
    {
        offsetMode_ = UNIFORM;
        offset_ = point(dict.lookup("offset"));
    }
    else if (dict.found("offsets"))
    {
        offsetMode_ = NONUNIFORM;
1088
1089
        //offsets_ = pointField(dict.lookup("offsets"));
        offsets_ = readListOrField("offsets", dict, patch_.size());
1090
    }
1091
    else if (mode_ != NEARESTPATCHFACE && mode_ != NEARESTPATCHFACEAMI)
1092
    {
1093
        FatalIOErrorInFunction
1094
1095
1096
        (
            dict
        )   << "Please supply the offsetMode as one of "
Mark OLESEN's avatar
Mark OLESEN committed
1097
            << offsetModeNames_
1098
1099
1100
1101
1102
            << exit(FatalIOError);
    }
}


1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
    const sampleMode mode,
    const dictionary& dict
)
:
    patch_(pp),
    sampleRegion_(dict.lookupOrDefault<word>("sampleRegion", "")),
    mode_(mode),
    samplePatch_(dict.lookupOrDefault<word>("samplePatch", "")),
    coupleGroup_(dict), //dict.lookupOrDefault<word>("coupleGroup", "")),
    offsetMode_(UNIFORM),
Henry Weller's avatar
Henry Weller committed
1116
    offset_(Zero),
1117
1118
1119
    offsets_(0),
    distance_(0.0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
1120
1121
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1122
    AMIReverse_(dict.lookupOrDefault<bool>("flipNormals", false)),
1123
    surfPtr_(nullptr),
1124
1125
1126
1127
    surfDict_(dict.subOrEmptyDict("surface"))
{
    if (mode != NEARESTPATCHFACE && mode != NEARESTPATCHFACEAMI)
    {
1128
        FatalIOErrorInFunction
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
        (
            dict
        )   << "Construct from sampleMode and dictionary only applicable for "
            << " collocated patches in modes "
            << sampleModeNames_[NEARESTPATCHFACE] << ','
            << sampleModeNames_[NEARESTPATCHFACEAMI]
            << exit(FatalIOError);
    }


    if (!coupleGroup_.valid())
    {
        if (sampleRegion_.empty())
        {
            // If no coupleGroup and no sampleRegion assume local region
            sampleRegion_ = patch_.boundaryMesh().mesh().name();
            sameRegion_ = true;
        }
    }
}


1151
1152
1153
Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
1154
    const mappedPatchBase& mpb
1155
1156
1157
)
:
    patch_(pp),
1158
1159
1160
    sampleRegion_(mpb.sampleRegion_),
    mode_(mpb.mode_),
    samplePatch_(mpb.samplePatch_),
1161
    coupleGroup_(mpb.coupleGroup_),
1162
1163
1164
1165
1166
    offsetMode_(mpb.offsetMode_),
    offset_(mpb.offset_),
    offsets_(mpb.offsets_),
    distance_(mpb.distance_),
    sameRegion_(mpb.sameRegion_),
1167
1168
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1169
    AMIReverse_(mpb.AMIReverse_),
1170
    surfPtr_(nullptr),
1171
    surfDict_(mpb.surfDict_)
1172
1173
1174
1175
1176
1177
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
1178
    const mappedPatchBase& mpb,
1179
1180
1181
1182
    const labelUList& mapAddressing
)
:
    patch_(pp),
1183
1184
1185
    sampleRegion_(mpb.sampleRegion_),
    mode_(mpb.mode_),
    samplePatch_(mpb.samplePatch_),
1186
    coupleGroup_(mpb.coupleGroup_),
1187
1188
    offsetMode_(mpb.offsetMode_),
    offset_(mpb.offset_),
1189
1190
1191
1192
1193
1194
    offsets_
    (
        offsetMode_ == NONUNIFORM
      ? vectorField(mpb.offsets_, mapAddressing)
      : vectorField(0)
    ),
1195
1196
    distance_(mpb.distance_),
    sameRegion_(mpb.sameRegion_),
1197
1198
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1199
    AMIReverse_(mpb.AMIReverse_),
1200
    surfPtr_(nullptr),
1201
    surfDict_(mpb.surfDict_)
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
{}


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

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


void Foam::mappedPatchBase::clearOut()
{
    mapPtr_.clear();
1216
1217
    AMIPtr_.clear();
    surfPtr_.clear();
1218
1219
1220
1221
1222
1223
1224
1225
1226
}


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

const Foam::polyMesh& Foam::mappedPatchBase::sampleMesh() const
{
    return patch_.boundaryMesh().mesh().time().lookupObject<polyMesh>
    (
1227
        sampleRegion()
1228
1229
1230
1231
1232
1233
1234
1235
    );
}


const Foam::polyPatch& Foam::mappedPatchBase::samplePolyPatch() const
{
    const polyMesh& nbrMesh = sampleMesh();

1236
    const label patchi = nbrMesh.boundaryMesh().findPatchID(samplePatch());
1237

1238
    if (patchi == -1)
1239
    {
1240
        FatalErrorInFunction
1241
            << "Cannot find patch " << samplePatch()
1242
1243
1244
1245
1246
            << " in region " << sampleRegion_ << endl
            << "Valid patches are " << nbrMesh.boundaryMesh().names()
            << exit(FatalError);
    }

1247
    return nbrMesh.boundaryMesh()[patchi];
1248
1249
1250
}


1251
1252
1253
1254
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints
(
    const pointField& fc
) const
1255
{
1256
    tmp<pointField> tfld(new pointField(fc));
1257
    pointField& fld = tfld.ref();
1258

1259
    switch (offsetMode_)
1260
1261
1262
1263
    {
        case UNIFORM:
        {
            fld += offset_;
1264
            break;
1265
1266
1267
1268
        }
        case NONUNIFORM:
        {
            fld += offsets_;
1269
            break;
1270
1271
1272
1273
1274
1275
1276
1277
        }
        case NORMAL:
        {
            // Get outwards pointing normal
            vectorField n(patch_.faceAreas());
            n /= mag(n);

            fld += distance_*n;
1278
            break;
1279
1280
1281
1282
1283
1284
1285
        }
    }

    return tfld;
}


1286
1287
1288
1289
1290
1291
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints() const
{
    return samplePoints(facePoints(patch_));
}


1292
1293
1294
Foam::pointIndexHit Foam::mappedPatchBase::facePoint
(
    const polyMesh& mesh,
1295
    const label facei,
1296
    const polyMesh::cellDecomposition decompMode
1297
1298
)
{
1299
    const point& fc = mesh.faceCentres()[facei];
1300
1301
1302

    switch (decompMode)
    {
1303
1304
        case polyMesh::FACE_PLANES:
        case polyMesh::FACE_CENTRE_TRIS:
1305
1306
1307
        {
            // For both decompositions the face centre is guaranteed to be
            // on the face
1308
            return pointIndexHit(true, fc, facei);
1309
1310
1311
        }
        break;

1312
1313
        case polyMesh::FACE_DIAG_TRIS:
        case polyMesh::CELL_TETS:
1314
        {
1315
            // Find the intersection of a ray from face centre to cell centre
1316
1317
1318
1319
            // Find intersection of (face-centre-decomposition) centre to
            // cell-centre with face-diagonal-decomposition triangles.

            const pointField& p = mesh.points();
1320
            const face& f = mesh.faces()[facei];
1321
1322
1323
1324
1325
1326
1327

            if (f.size() <= 3)
            {
                // Return centre of triangle.
                return pointIndexHit(true, fc, 0);
            }

1328
1329
            label celli = mesh.faceOwner()[facei];
            const point& cc = mesh.cellCentres()[celli];
1330
1331
            vector d = fc-cc;

1332
            const label fp0 = mesh.tetBasePtIs()[facei];
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
            const point& basePoint = p[f[fp0]];

            label fp = f.fcIndex(fp0);
            for (label i = 2; i < f.size(); i++)
            {
                const point& thisPoint = p[f[fp]];
                label nextFp = f.fcIndex(fp);
                const point& nextPoint = p[f[nextFp]];

                const triPointRef tri(basePoint, thisPoint, nextPoint);
                pointHit hitInfo = tri.intersection
                (
                    cc,
                    d,
                    intersection::HALF_RAY
                );

                if (hitInfo.hit() && hitInfo.distance() > 0)
                {
                    return pointIndexHit(true, hitInfo.hitPoint(), i-2);
                }

                fp = nextFp;
            }

            // Fall-back
            return pointIndexHit(false, fc, -1);
        }
        break;

        default:
        {
1365
            FatalErrorInFunction
1366
1367
1368
1369
1370
1371
1372
                << "problem" << abort(FatalError);
            return pointIndexHit();
        }
    }
}


1373
1374
void Foam::mappedPatchBase::write(Ostream& os) const
{
1375
    os.writeEntry("sampleMode", sampleModeNames_[mode_]);
1376
1377
    if (!sampleRegion_.empty())
    {
1378
        os.writeEntry("sampleRegion", sampleRegion_);
1379
1380
1381
    }
    if (!samplePatch_.empty())
    {
1382
        os.writeEntry("samplePatch", samplePatch_);
1383
    }
1384
    coupleGroup_.write(os);
1385

1386
1387
1388
1389
1390
1391
    if
    (
        offsetMode_ == UNIFORM
     && offset_ == vector::zero
     && (mode_ == NEARESTPATCHFACE || mode_ == NEARESTPATCHFACEAMI)
    )
1392
    {
1393
        // Collocated mode. No need to write offset data
1394
    }
1395
    else
1396
    {
1397
        os.writeEntry("offsetMode", offsetModeNames_[offsetMode_]);
1398
1399

        switch (offsetMode_)
1400
        {
1401
1402
            case UNIFORM:
            {
1403
                os.writeEntry("offset", offset_);
1404
1405
1406
1407
1408
1409
1410
1411
1412
                break;
            }
            case NONUNIFORM:
            {
                offsets_.writeEntry("offsets", os);
                break;
            }
            case NORMAL:
            {