mappedPatchBase.C 39.4 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
OpenFOAM bot's avatar
OpenFOAM bot committed
6
7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
9
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2015-2018 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"
48
49
50
51
52
53
54
55
56

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

namespace Foam
{
    defineTypeNameAndDebug(mappedPatchBase, 0);
}


57
58
59
60
61
const Foam::Enum
<
    Foam::mappedPatchBase::sampleMode
>
Foam::mappedPatchBase::sampleModeNames_
Mark OLESEN's avatar
Mark OLESEN committed
62
({
63
64
65
66
67
68
    { 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
69
});
70
71
72
73
74
75
76


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


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

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

    // Force construction of min-tet decomp
    (void)mesh.tetBasePtIs();
95
96

    // Initialise to face-centre
97
    tmp<pointField> tfacePoints(new pointField(patch_.size()));
98
    pointField& facePoints = tfacePoints.ref();
99

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

    return tfacePoints;
}


114
115
void Foam::mappedPatchBase::collectSamples
(
116
    const pointField& facePoints,
117
118
119
120
121
122
123
    pointField& samples,
    labelList& patchFaceProcs,
    labelList& patchFaces,
    pointField& patchFc
) const
{
    // Collect all sample points and the faces they come from.
124
125
126
127
128
129
130
131
132
133
134
135
    {
        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>()
        );
    }
136

137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
    {
        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>
158
159
160
        (
            globalFaces,
            accessOp<labelList>()
161
162
163
        );
    }

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

        patchFaceProcs.setSize(patchFaces.size());

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


// Find the processor/cell containing the samples. Does not account
// for samples being found in two processors.
void Foam::mappedPatchBase::findSamples
(
188
    const sampleMode mode,
189
190
191
192
193
194
195
196
197
198
199
200
    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());

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

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

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

219
                label celli = tree.findInside(sample);
220

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

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

243
244
245
246
        case NEARESTONLYCELL:
        {
            if (samplePatch_.size() && samplePatch_ != "none")
            {
247
248
                FatalErrorInFunction
                    << "No need to supply a patch name when in "
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
                    << 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;
        }

270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
        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
287
                const labelList patchFaces(identity(pp.size(), pp.start()));
288
289
290
291
292
293

                treeBoundBox patchBb
                (
                    treeBoundBox(pp.points(), pp.meshPoints()).extend
                    (
                        rndGen,
294
                        1e-4
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
341
342
343
                    )
                );
                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;
        }

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
412
413
414
        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;
        }

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

424
425
426
            //- Note: face-diagonal decomposition
            const meshSearchMeshObject& meshSearchEngine =
                meshSearchMeshObject::New(mesh);
427
428
429
430
431

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

432
                label facei = meshSearchEngine.findNearestFace(sample);
433

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

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

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

        default:
        {
464
            FatalErrorInFunction
465
466
467
468
469
                << "problem." << abort(FatalError);
        }
    }


470
    // Find nearest. Combine on master.
471
472
473
    Pstream::listCombineGather(nearest, nearestEqOp());
    Pstream::listCombineScatter(nearest);

474

475
476
    if (debug)
    {
477
478
479
        InfoInFunction
            << "mesh " << sampleRegion() << " : " << endl;

480
481
        forAll(nearest, sampleI)
        {
482
            label proci = nearest[sampleI].second().second();
483
484
485
            label localI = nearest[sampleI].first().index();

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

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

    forAll(nearest, sampleI)
    {
500
501
502
503
504
505
506
507
508
509
510
511
        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();
        }
512
513
514
515
516
517
    }
}


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

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

    // Get offsetted points
530
    const pointField offsettedPoints(samplePoints(patchPoints()));
531

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

    // Check offset
542
543
544
545
    vectorField d(offsettedPoints-patchPoints());
    bool coincident = (gAverage(mag(d)) <= ROOTVSMALL);

    if (sampleMyself && coincident)
546
    {
547
548
        WarningInFunction
            << "Invalid offset " << d << endl
549
550
551
552
553
554
555
            << "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
556
            << "sampleRegion_:" << sampleRegion() << endl
557
            << "mode_:" << sampleModeNames_[mode_] << endl
558
            << "samplePatch_:" << samplePatch() << endl
559
560
561
562
563
564
565
566
            << "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;
567
568
569
570
571
572
573
574
575
    collectSamples
    (
        patchPoints,
        samples,
        patchFaceProcs,
        patchFaces,
        patchFc
    );

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

582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
    // 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)
            {
600
601
                WarningInFunction
                    << "Did not find " << nNotFound
602
603
                    << " out of " << sampleProcs.size() << " total samples."
                    << " Sampling these on owner cell centre instead." << endl
604
                    << "On patch " << patch_.name()
605
                    << " on region " << sampleRegion()
606
                    << " in mode " << sampleModeNames_[mode_] << endl
607
608
                    << "with offset mode " << offsetModeNames_[offsetMode_]
                    << ". Suppressing further warnings from " << type() << endl;
609
610
611
612

                hasWarned = true;
            }

613
614
615
            // Collect the samples that cannot be found
            DynamicList<label> subMap;
            DynamicField<point> subSamples;
616
617
618
619
620

            forAll(sampleProcs, sampleI)
            {
                if (sampleProcs[sampleI] == -1)
                {
621
622
                    subMap.append(sampleI);
                    subSamples.append(samples[sampleI]);
623
624
625
                }
            }

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

            // Insert
640
641
            labelUIndList(sampleProcs, subMap) = subSampleProcs;
            labelUIndList(sampleIndices, subMap) = subSampleIndices;
642
            UIndirectList<point>(sampleLocations, subMap) = subSampleLocations;
643
644
645
        }
    }

646
647
648
649
650
651
652
653
654
655
656
657
658
    // 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
659
    //        << "Found data in region " << sampleRegion()
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
    //        << " 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"
677
678
            << " sampled cell/faceCentres/points to file " << str.name()
            << endl;
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700

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

701
    forAll(subMap, proci)
702
    {
703
704
        subMap[proci] = labelUIndList(sampleIndices, subMap[proci]);
        constructMap[proci] = labelUIndList(patchFaces, constructMap[proci]);
705
706
707

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

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

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

            forAll(map, i)
            {
729
                label facei = map[i];
730

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


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

    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())
    {
796
        FatalErrorInFunction
797
798
799
800
801
            << "AMI already calculated" << exit(FatalError);
    }

    AMIPtr_.clear();

802

803
    const polyPatch& nbr = samplePolyPatch();
804

805
806
    // Transform neighbour patch to local system
    pointField nbrPoints(samplePoints(nbr.localPoints()));
807

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


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

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

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


849
850
851
852
853
854
855
856
857
858
859
// 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());
860
    pointField& fld = tfld.ref();
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880

    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)
                {
881
882
                    FatalIOErrorInFunction(dict)
                        << "size " << fld.size()
883
884
885
886
887
888
                        << " is not equal to the given value of " << size
                        << exit(FatalIOError);
                }
            }
            else
            {
889
                FatalIOErrorInFunction(dict)
890
                    << "Expected keyword 'uniform' or 'nonuniform', found "
891
892
893
894
                    << firstToken.wordToken()
                    << exit(FatalIOError);
            }
        }
895
        else if (is.version() == IOstream::versionNumber(2,0))
896
        {
897
            IOWarningInFunction(dict)
898
                << "Expected keyword 'uniform' or 'nonuniform', "
899
900
                   "assuming List format for backwards compatibility."
                   "Foam version 2.0." << endl;
901

902
903
            is.putBack(firstToken);
            is >> static_cast<List<vector>&>(fld);
904
905
906
907
908
909
        }
    }
    return tfld;
}


910
911
912
913
914
915
916
917
918
919
// * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //

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


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),
948
    coupleGroup_(),
949
    offsetMode_(NONUNIFORM),
Henry Weller's avatar
Henry Weller committed
950
    offset_(Zero),
951
952
953
    offsets_(offsets),
    distance_(0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
954
955
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
956
    AMIReverse_(false),
957
    surfPtr_(nullptr),
958
    surfDict_(fileName("surface"))
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
{}


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),
975
    coupleGroup_(),
976
977
978
979
980
    offsetMode_(UNIFORM),
    offset_(offset),
    offsets_(0),
    distance_(0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
981
982
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
983
    AMIReverse_(false),
984
    surfPtr_(nullptr),
985
    surfDict_(fileName("surface"))
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
{}


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),
1002
    coupleGroup_(),
1003
    offsetMode_(NORMAL),
Henry Weller's avatar
Henry Weller committed
1004
    offset_(Zero),
1005
1006
1007
    offsets_(0),
    distance_(distance),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
1008
1009
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1010
    AMIReverse_(false),
1011
    surfPtr_(nullptr),
1012
    surfDict_(fileName("surface"))
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
{}


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

Mark OLESEN's avatar
Mark OLESEN committed
1048
    if (offsetModeNames_.readIfPresent("offsetMode", dict, offsetMode_))
1049
    {
Mark OLESEN's avatar
Mark OLESEN committed
1050
        switch (offsetMode_)
1051
1052
1053
        {
            case UNIFORM:
            {
1054
                dict.readEntry("offset", offset_);
1055
1056
1057
1058
1059
            }
            break;

            case NONUNIFORM:
            {
1060
1061
                //offsets_ = pointField(dict.lookup("offsets"));
                offsets_ = readListOrField("offsets", dict, patch_.size());
1062
1063
1064
1065
1066
            }
            break;

            case NORMAL:
            {
1067
                dict.readEntry("distance", distance_);
1068
1069
1070
1071
            }
            break;
        }
    }
1072
    else if (dict.readIfPresent("offset", offset_))
1073
1074
1075
1076
1077
1078
    {
        offsetMode_ = UNIFORM;
    }
    else if (dict.found("offsets"))
    {
        offsetMode_ = NONUNIFORM;
1079
1080
        //offsets_ = pointField(dict.lookup("offsets"));
        offsets_ = readListOrField("offsets", dict, patch_.size());
1081
    }
1082
    else if (mode_ != NEARESTPATCHFACE && mode_ != NEARESTPATCHFACEAMI)
1083
    {
1084
1085
        FatalIOErrorInFunction(dict)
            << "Please supply the offsetMode as one of "
Mark OLESEN's avatar
Mark OLESEN committed
1086
            << offsetModeNames_
1087
1088
1089
1090
1091
            << exit(FatalIOError);
    }
}


1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
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
1105
    offset_(Zero),
1106
1107
1108
    offsets_(0),
    distance_(0.0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
1109
1110
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1111
    AMIReverse_(dict.lookupOrDefault("flipNormals", false)),
1112
    surfPtr_(nullptr),
1113
1114
1115
1116
    surfDict_(dict.subOrEmptyDict("surface"))
{
    if (mode != NEARESTPATCHFACE && mode != NEARESTPATCHFACEAMI)
    {
1117
1118
        FatalIOErrorInFunction(dict)
            << "Construct from sampleMode and dictionary only applicable for "
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
            << " 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;
        }
    }
}


1138
1139
1140
Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
1141
    const mappedPatchBase& mpb
1142
1143
1144
)
:
    patch_(pp),
1145
1146
1147
    sampleRegion_(mpb.sampleRegion_),
    mode_(mpb.mode_),
    samplePatch_(mpb.samplePatch_),
1148
    coupleGroup_(mpb.coupleGroup_),
1149
1150
1151
1152
1153
    offsetMode_(mpb.offsetMode_),
    offset_(mpb.offset_),
    offsets_(mpb.offsets_),
    distance_(mpb.distance_),
    sameRegion_(mpb.sameRegion_),
1154
1155
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1156
    AMIReverse_(mpb.AMIReverse_),
1157
    surfPtr_(nullptr),
1158
    surfDict_(mpb.surfDict_)
1159
1160
1161
1162
1163
1164
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
1165
    const mappedPatchBase& mpb,
1166
1167
1168
1169
    const labelUList& mapAddressing
)
:
    patch_(pp),
1170
1171
1172
    sampleRegion_(mpb.sampleRegion_),
    mode_(mpb.mode_),
    samplePatch_(mpb.samplePatch_),
1173
    coupleGroup_(mpb.coupleGroup_),
1174
1175
    offsetMode_(mpb.offsetMode_),
    offset_(mpb.offset_),
1176
1177
1178
1179
1180
1181
    offsets_
    (
        offsetMode_ == NONUNIFORM
      ? vectorField(mpb.offsets_, mapAddressing)
      : vectorField(0)
    ),
1182
1183
    distance_(mpb.distance_),
    sameRegion_(mpb.sameRegion_),
1184
1185
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1186
    AMIReverse_(mpb.AMIReverse_),
1187
    surfPtr_(nullptr),
1188
    surfDict_(mpb.surfDict_)
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
{}


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

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


void Foam::mappedPatchBase::clearOut()
{
    mapPtr_.clear();
1203
1204
    AMIPtr_.clear();
    surfPtr_.clear();
1205
1206
1207
1208
1209
1210
1211
}


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

const Foam::polyMesh& Foam::mappedPatchBase::sampleMesh() const
{
1212
1213
1214
    const polyMesh& thisMesh = patch_.boundaryMesh().mesh();

    return
1215
    (
1216
1217
1218
        sameRegion_
      ? thisMesh
      : thisMesh.time().lookupObject<polyMesh>(sampleRegion())
1219
1220
1221
1222
1223
1224
1225
1226
    );
}


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

1227
    const label patchi = nbrMesh.boundaryMesh().findPatchID(samplePatch());
1228

1229
    if (patchi == -1)
1230
    {
1231
        FatalErrorInFunction
1232
            << "Cannot find patch " << samplePatch()
1233
1234
1235
1236
1237
            << " in region " << sampleRegion_ << endl
            << "Valid patches are " << nbrMesh.boundaryMesh().names()
            << exit(FatalError);
    }

1238
    return nbrMesh.boundaryMesh()[patchi];
1239
1240
1241
}


1242
1243
1244
1245
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints
(
    const pointField& fc
) const
1246
{
1247
    tmp<pointField> tfld(new pointField(fc));
1248
    pointField& fld = tfld.ref();
1249

1250
    switch (offsetMode_)
1251
1252
1253
1254
    {
        case UNIFORM:
        {
            fld += offset_;
1255
            break;
1256
1257
1258
1259
        }
        case NONUNIFORM:
        {
            fld += offsets_;
1260
            break;
1261
1262
1263
1264
1265
1266
1267
1268
        }
        case NORMAL:
        {
            // Get outwards pointing normal
            vectorField n(patch_.faceAreas());
            n /= mag(n);

            fld += distance_*n;
1269
            break;
1270
1271
1272
1273
1274
1275
1276
        }
    }

    return tfld;
}


1277
1278
1279
1280
1281
1282
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints() const
{
    return samplePoints(facePoints(patch_));
}


1283
1284
1285
Foam::pointIndexHit Foam::mappedPatchBase::facePoint
(
    const polyMesh& mesh,
1286
    const label facei,
1287
    const polyMesh::cellDecomposition decompMode
1288
1289
)
{
1290
    const point& fc = mesh.faceCentres()[facei];
1291
1292
1293

    switch (decompMode)
    {
1294
1295
        case polyMesh::FACE_PLANES:
        case polyMesh::FACE_CENTRE_TRIS:
1296
1297
1298
        {
            // For both decompositions the face centre is guaranteed to be
            // on the face
1299
            return pointIndexHit(true, fc, facei);
1300
1301
1302
        }
        break;

1303
1304
        case polyMesh::FACE_DIAG_TRIS:
        case polyMesh::CELL_TETS:
1305
        {
1306
            // Find the intersection of a ray from face centre to cell centre
1307
1308
1309
1310
            // Find intersection of (face-centre-decomposition) centre to
            // cell-centre with face-diagonal-decomposition triangles.

            const pointField& p = mesh.points();
1311
            const face& f = mesh.faces()[facei];
1312
1313
1314
1315
1316
1317
1318

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

1319
1320
            label celli = mesh.faceOwner()[facei];
            const point& cc = mesh.cellCentres()[celli];
1321
1322
            vector d = fc-cc;

1323
            const label fp0 = mesh.tetBasePtIs()[facei];
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
            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:
        {
1356
            FatalErrorInFunction
1357
1358
1359
1360
1361
1362
1363
                << "problem" << abort(FatalError);
            return pointIndexHit();
        }
    }
}


1364
1365
void Foam::mappedPatchBase::write(Ostream& os) const
{
1366
    os.writeEntry("sampleMode", sampleModeNames_[mode_]);
1367
1368
    if (!sampleRegion_.empty())
    {
1369
        os.writeEntry("sampleRegion", sampleRegion_);
1370
1371
1372
    }
    if (!samplePatch_.empty())
    {
1373
        os.writeEntry("samplePatch", samplePatch_);
1374
    }