mappedPatchBase.C 39.6 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
284
285
286
287
288
289
290
        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
                const labelList patchFaces(identity(pp.size()) + pp.start());

                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
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
            UIndirectList<label>(sampleProcs, subMap) = subSampleProcs;
            UIndirectList<label>(sampleIndices, subMap) = subSampleIndices;
            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
        subMap[proci] = UIndirectList<label>
701
702
        (
            sampleIndices,
703
            subMap[proci]
704
        );
705
        constructMap[proci] = UIndirectList<label>
706
707
        (
            patchFaces,
708
            constructMap[proci]
709
710
711
712
        );

        //if (debug)
        //{
713
714
715
        //    Pout<< "To proc:" << proci << " sending values of cells/faces:"
        //        << subMap[proci] << endl;
        //    Pout<< "From proc:" << proci
716
        //        << " receiving values of patch faces:"
717
        //        << constructMap[proci] << endl;
718
719
720
721
722
723
724
725
726
727
        //}
    }

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

    if (debug)
    {
        // Check that all elements get a value.
        PackedBoolList used(patch_.size());
728
        forAll(constructMap, proci)
729
        {
730
            const labelList& map = constructMap[proci];
731
732
733

            forAll(map, i)
            {
734
                label facei = map[i];
735

736
                if (used[facei] == 0)
737
                {
738
                    used[facei] = 1;
739
740
741
                }
                else
                {
742
                    FatalErrorInFunction
743
                        << "On patch " << patch_.name()
744
                        << " patchface " << facei
745
746
747
748
749
                        << " is assigned to more than once."
                        << abort(FatalError);
                }
            }
        }
750
        forAll(used, facei)
751
        {
752
            if (used[facei] == 0)
753
            {
754
                FatalErrorInFunction
755
                    << "On patch " << patch_.name()
756
                    << " patchface " << facei
757
758
759
760
761
762
763
764
765
766
767
                    << " is never assigned to."
                    << abort(FatalError);
            }
        }
    }
}


const Foam::autoPtr<Foam::searchableSurface>& Foam::mappedPatchBase::surfPtr()
const
{
768
    const word surfType(surfDict_.lookupOrDefault<word>("type", "none"));
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
796
797
798
799
800

    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())
    {
801
        FatalErrorInFunction
802
803
804
805
806
            << "AMI already calculated" << exit(FatalError);
    }

    AMIPtr_.clear();

807

808
    const polyPatch& nbr = samplePolyPatch();
809

810
811
    // Transform neighbour patch to local system
    pointField nbrPoints(samplePoints(nbr.localPoints()));
812

813
814
815
    primitivePatch nbrPatch0
    (
        SubList<face>
816
        (
817
818
819
820
821
822
823
824
825
826
827
            nbr.localFaces(),
            nbr.size()
        ),
        nbrPoints
    );


    if (debug)
    {
        OFstream os(patch_.name() + "_neighbourPatch-org.obj");
        meshTools::writeOBJ(os, samplePolyPatch().localFaces(), nbrPoints);
828
829
830
831
832
833
834

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

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


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

    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)
                {
886
                    FatalIOErrorInFunction
887
888
889
890
891
892
893
894
895
                    (
                        dict
                    )   << "size " << fld.size()
                        << " is not equal to the given value of " << size
                        << exit(FatalIOError);
                }
            }
            else
            {
896
                FatalIOErrorInFunction
897
898
899
900
901
902
903
904
905
906
907
                (
                    dict
                )   << "expected keyword 'uniform' or 'nonuniform', found "
                    << firstToken.wordToken()
                    << exit(FatalIOError);
            }
        }
        else
        {
            if (is.version() == 2.0)
            {
908
                IOWarningInFunction
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
                (
                    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;
}


924
925
926
927
928
929
930
931
932
933
// * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //

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


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),
962
    coupleGroup_(),
963
    offsetMode_(NONUNIFORM),
Henry Weller's avatar
Henry Weller committed
964
    offset_(Zero),
965
966
967
    offsets_(offsets),
    distance_(0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
968
969
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
970
    AMIReverse_(false),
971
    surfPtr_(nullptr),
972
    surfDict_(fileName("surface"))
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
{}


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),
989
    coupleGroup_(),
990
991
992
993
994
    offsetMode_(UNIFORM),
    offset_(offset),
    offsets_(0),
    distance_(0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
995
996
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
997
    AMIReverse_(false),
998
    surfPtr_(nullptr),
999
    surfDict_(fileName("surface"))
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
{}


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),
1016
    coupleGroup_(),
1017
    offsetMode_(NORMAL),
Henry Weller's avatar
Henry Weller committed
1018
    offset_(Zero),
1019
1020
1021
    offsets_(0),
    distance_(distance),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
1022
1023
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1024
    AMIReverse_(false),
1025
    surfPtr_(nullptr),
1026
    surfDict_(fileName("surface"))
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
{}


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

1062
1063
    if (dict.found("offsetMode"))
    {
1064
        offsetMode_ = offsetModeNames_.lookup("offsetMode", dict);
1065

Mark OLESEN's avatar
Mark OLESEN committed
1066
        switch (offsetMode_)
1067
1068
1069
1070
1071
1072
1073
1074
1075
        {
            case UNIFORM:
            {
                offset_ = point(dict.lookup("offset"));
            }
            break;

            case NONUNIFORM:
            {
1076
1077
                //offsets_ = pointField(dict.lookup("offsets"));
                offsets_ = readListOrField("offsets", dict, patch_.size());
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
            }
            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;
1096
1097
        //offsets_ = pointField(dict.lookup("offsets"));
        offsets_ = readListOrField("offsets", dict, patch_.size());
1098
    }
1099
    else if (mode_ != NEARESTPATCHFACE && mode_ != NEARESTPATCHFACEAMI)
1100
    {
1101
        FatalIOErrorInFunction
1102
1103
1104
        (
            dict
        )   << "Please supply the offsetMode as one of "
Mark OLESEN's avatar
Mark OLESEN committed
1105
            << offsetModeNames_
1106
1107
1108
1109
1110
            << exit(FatalIOError);
    }
}


1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
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
1124
    offset_(Zero),
1125
1126
1127
    offsets_(0),
    distance_(0.0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
1128
1129
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1130
    AMIReverse_(dict.lookupOrDefault<bool>("flipNormals", false)),
1131
    surfPtr_(nullptr),
1132
1133
1134
1135
    surfDict_(dict.subOrEmptyDict("surface"))
{
    if (mode != NEARESTPATCHFACE && mode != NEARESTPATCHFACEAMI)
    {
1136
        FatalIOErrorInFunction
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
        (
            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;
        }
    }
}


1159
1160
1161
Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
1162
    const mappedPatchBase& mpb
1163
1164
1165
)
:
    patch_(pp),
1166
1167
1168
    sampleRegion_(mpb.sampleRegion_),
    mode_(mpb.mode_),
    samplePatch_(mpb.samplePatch_),
1169
    coupleGroup_(mpb.coupleGroup_),
1170
1171
1172
1173
1174
    offsetMode_(mpb.offsetMode_),
    offset_(mpb.offset_),
    offsets_(mpb.offsets_),
    distance_(mpb.distance_),
    sameRegion_(mpb.sameRegion_),
1175
1176
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1177
    AMIReverse_(mpb.AMIReverse_),
1178
    surfPtr_(nullptr),
1179
    surfDict_(mpb.surfDict_)
1180
1181
1182
1183
1184
1185
{}


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


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

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


void Foam::mappedPatchBase::clearOut()
{
    mapPtr_.clear();
1224
1225
    AMIPtr_.clear();
    surfPtr_.clear();
1226
1227
1228
1229
1230
1231
1232
1233
1234
}


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

const Foam::polyMesh& Foam::mappedPatchBase::sampleMesh() const
{
    return patch_.boundaryMesh().mesh().time().lookupObject<polyMesh>
    (
1235
        sampleRegion()
1236
1237
1238
1239
1240
1241
1242
1243
    );
}


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

1244
    const label patchi = nbrMesh.boundaryMesh().findPatchID(samplePatch());
1245

1246
    if (patchi == -1)
1247
    {
1248
        FatalErrorInFunction
1249
            << "Cannot find patch " << samplePatch()
1250
1251
1252
1253
1254
            << " in region " << sampleRegion_ << endl
            << "Valid patches are " << nbrMesh.boundaryMesh().names()
            << exit(FatalError);
    }

1255
    return nbrMesh.boundaryMesh()[patchi];
1256
1257
1258
}


1259
1260
1261
1262
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints
(
    const pointField& fc
) const
1263
{
1264
    tmp<pointField> tfld(new pointField(fc));
1265
    pointField& fld = tfld.ref();
1266

1267
    switch (offsetMode_)
1268
1269
1270
1271
    {
        case UNIFORM:
        {
            fld += offset_;
1272
            break;
1273
1274
1275
1276
        }
        case NONUNIFORM:
        {
            fld += offsets_;
1277
            break;
1278
1279
1280
1281
1282
1283
1284
1285
        }
        case NORMAL:
        {
            // Get outwards pointing normal
            vectorField n(patch_.faceAreas());
            n /= mag(n);

            fld += distance_*n;
1286
            break;
1287
1288
1289
1290
1291
1292
1293
        }
    }

    return tfld;
}


1294
1295
1296
1297
1298
1299
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints() const
{
    return samplePoints(facePoints(patch_));
}


1300
1301
1302
Foam::pointIndexHit Foam::mappedPatchBase::facePoint
(
    const polyMesh& mesh,
1303
    const label facei,
1304
    const polyMesh::cellDecomposition decompMode
1305
1306
)
{
1307
    const point& fc = mesh.faceCentres()[facei];
1308
1309
1310

    switch (decompMode)
    {
1311
1312
        case polyMesh::FACE_PLANES:
        case polyMesh::FACE_CENTRE_TRIS:
1313
1314
1315
        {
            // For both decompositions the face centre is guaranteed to be
            // on the face
1316
            return pointIndexHit(true, fc, facei);
1317
1318
1319
        }
        break;

1320
1321
        case polyMesh::FACE_DIAG_TRIS:
        case polyMesh::CELL_TETS:
1322
        {
1323
            // Find the intersection of a ray from face centre to cell centre
1324
1325
1326
1327
            // Find intersection of (face-centre-decomposition) centre to
            // cell-centre with face-diagonal-decomposition triangles.

            const pointField& p = mesh.points();
1328
            const face& f = mesh.faces()[facei];
1329
1330
1331
1332
1333
1334
1335

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

1336
1337
            label celli = mesh.faceOwner()[facei];
            const point& cc = mesh.cellCentres()[celli];
1338
1339
            vector d = fc-cc;

1340
            const label fp0 = mesh.tetBasePtIs()[facei];
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
            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:
        {
1373
            FatalErrorInFunction
1374
1375
1376
1377
1378
1379
1380
                << "problem" << abort(FatalError);
            return pointIndexHit();
        }
    }
}


1381
1382
void Foam::mappedPatchBase::write(Ostream& os) const
{
1383
    os.writeEntry("sampleMode", sampleModeNames_[mode_]);
1384
1385
    if (!sampleRegion_.empty())
    {
1386
        os.writeEntry("sampleRegion", sampleRegion_);
1387
1388
1389
    }
    if (!samplePatch_.empty())
    {
1390
        os.writeEntry("samplePatch", samplePatch_);
1391
    }
1392
    coupleGroup_.write(os);
1393

1394
1395
1396
1397
1398
1399
    if
    (
        offsetMode_ == UNIFORM
     && offset_ == vector::zero
     && (mode_ == NEARESTPATCHFACE || mode_ == NEARESTPATCHFACEAMI)
    )
1400
    {
1401
        // Collocated mode. No need to write offset data
1402
    }
1403
    else
1404
    {
1405
        os.writeEntry("offsetMode", offsetModeNames_[offsetMode_]);
1406
1407

        switch (offsetMode_)
1408
        {
1409
1410
            case UNIFORM:
            {
1411
                os.writeEntry("offset", offset_);