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


const Foam::Enum
<
    Foam::mappedPatchBase::offsetMode
>
Foam::mappedPatchBase::offsetModeNames_
Mark OLESEN's avatar
Mark OLESEN committed
74
({
75
76
77
    { offsetMode::UNIFORM, "uniform" },
    { offsetMode::NONUNIFORM, "nonuniform" },
    { offsetMode::NORMAL, "normal" },
Mark OLESEN's avatar
Mark OLESEN committed
78
});
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
        //}
    }

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

    if (debug)
    {
        // Check that all elements get a value.
719
        bitSet 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.test(facei))
729
                {
730
                    FatalErrorInFunction
731
                        << "On patch " << patch_.name()
732
                        << " patchface " << facei
733
734
735
                        << " is assigned to more than once."
                        << abort(FatalError);
                }
736
737
738
739
                else
                {
                    used.set(facei);
                }
740
741
            }
        }
742
        forAll(used, facei)
743
        {
744
            if (!used.test(facei))
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
879
                    FatalIOErrorInFunction(dict)
                        << "size " << fld.size()
880
881
882
883
884
885
                        << " is not equal to the given value of " << size
                        << exit(FatalIOError);
                }
            }
            else
            {
886
887
                FatalIOErrorInFunction(dict)
                    << "expected keyword 'uniform' or 'nonuniform', found "
888
889
890
891
                    << firstToken.wordToken()
                    << exit(FatalIOError);
            }
        }
892
        else if (is.version() == IOstream::versionNumber(2,0))
893
        {
894
895
896
897
            IOWarningInFunction(dict)
                << "expected keyword 'uniform' or 'nonuniform', "
                   "assuming List format for backwards compatibility."
                   "Foam version 2.0." << endl;
898

899
900
            is.putBack(firstToken);
            is >> static_cast<List<vector>&>(fld);
901
902
903
904
905
906
        }
    }
    return tfld;
}


907
908
909
910
911
912
913
914
915
916
// * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //

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


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


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


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


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
    const dictionary& dict
)
:
    patch_(pp),
1020
    sampleRegion_(dict.lookupOrDefault<word>("sampleRegion", "")),
Mark OLESEN's avatar
Mark OLESEN committed
1021
    mode_(sampleModeNames_.get("sampleMode", dict)),
1022
1023
    samplePatch_(dict.lookupOrDefault<word>("samplePatch", "")),
    coupleGroup_(dict),
1024
    offsetMode_(UNIFORM),
Henry Weller's avatar
Henry Weller committed
1025
    offset_(Zero),
1026
1027
1028
    offsets_(0),
    distance_(0.0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
1029
1030
    mapPtr_(nullptr),
    AMIPtr_(nullptr),
1031
    AMIReverse_(dict.lookupOrDefault("flipNormals", false)),
1032
    surfPtr_(nullptr),
1033
1034
    surfDict_(dict.subOrEmptyDict("surface"))
{
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
    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
1045
    if (offsetModeNames_.readIfPresent("offsetMode", dict, offsetMode_))
1046
    {
Mark OLESEN's avatar
Mark OLESEN committed
1047
        switch (offsetMode_)
1048
1049
1050
        {
            case UNIFORM:
            {
1051
                dict.readEntry("offset", offset_);
1052
1053
1054
1055
1056
            }
            break;

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

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


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


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


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


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

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


void Foam::mappedPatchBase::clearOut()
{
    mapPtr_.clear();
1200
1201
    AMIPtr_.clear();
    surfPtr_.clear();
1202
1203
1204
1205
1206
1207
1208
}


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

const Foam::polyMesh& Foam::mappedPatchBase::sampleMesh() const
{
1209
1210
1211
    const polyMesh& thisMesh = patch_.boundaryMesh().mesh();

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


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

1224
    const label patchi = nbrMesh.boundaryMesh().findPatchID(samplePatch());
1225

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

1235
    return nbrMesh.boundaryMesh()[patchi];
1236
1237
1238
}


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

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

            fld += distance_*n;
1266
            break;
1267
1268
1269
1270
1271
1272
1273
        }
    }

    return tfld;
}


1274
1275
1276
1277
1278
1279
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints() const
{
    return samplePoints(facePoints(patch_));
}


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

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

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

            const pointField& p = mesh.points();
1308
            const face& f = mesh.faces()[facei];
1309
1310
1311
1312
1313
1314
1315

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

1316
1317
            label celli = mesh.faceOwner()[facei];
            const point& cc = mesh.cellCentres()[celli];
1318
1319
            vector d = fc-cc;

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


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

1374
1375
1376
1377
1378
1379
    if
    (
        offsetMode_ == UNIFORM
     && offset_ == vector::zero
     && (mode_ == NEARESTPATCHFACE || mode_ == NEARESTPATCHFACEAMI)
    )
andy's avatar