mappedPatchBase.C 42 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2015 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
54
55

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

namespace Foam
{
    defineTypeNameAndDebug(mappedPatchBase, 0);

    template<>
    const char* Foam::NamedEnum
    <
        Foam::mappedPatchBase::sampleMode,
56
        6
57
58
59
60
    >::names[] =
    {
        "nearestCell",
        "nearestPatchFace",
61
        "nearestPatchFaceAMI",
62
        "nearestPatchPoint",
63
64
        "nearestFace",
        "nearestOnlyCell"
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
    };

    template<>
    const char* Foam::NamedEnum
    <
        Foam::mappedPatchBase::offsetMode,
        3
    >::names[] =
    {
        "uniform",
        "nonuniform",
        "normal"
    };
}


81
const Foam::NamedEnum<Foam::mappedPatchBase::sampleMode, 6>
82
83
84
85
86
87
88
89
    Foam::mappedPatchBase::sampleModeNames_;

const Foam::NamedEnum<Foam::mappedPatchBase::offsetMode, 3>
    Foam::mappedPatchBase::offsetModeNames_;


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

90
91
92
93
94
95
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::facePoints
(
    const polyPatch& pp
) const
{
    const polyMesh& mesh = pp.boundaryMesh().mesh();
96
97
98

    // Force construction of min-tet decomp
    (void)mesh.tetBasePtIs();
99
100

    // Initialise to face-centre
101
    tmp<pointField> tfacePoints(new pointField(patch_.size()));
102
103
104
105
    pointField& facePoints = tfacePoints();

    forAll(pp, faceI)
    {
106
107
108
109
        facePoints[faceI] = facePoint
        (
            mesh,
            pp.start()+faceI,
110
            polyMesh::FACE_DIAG_TRIS
111
        ).rawPoint();
112
113
114
115
116
117
    }

    return tfacePoints;
}


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

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

168
    {
169
170
171
172
173
174
175
176
177
        labelList nPerProc(Pstream::nProcs());
        nPerProc[Pstream::myProcNo()] = patch_.size();
        Pstream::gatherList(nPerProc);
        Pstream::scatterList(nPerProc);

        patchFaceProcs.setSize(patchFaces.size());

        label sampleI = 0;
        forAll(nPerProc, procI)
178
        {
179
180
181
182
            for (label i = 0; i < nPerProc[procI]; i++)
            {
                patchFaceProcs[sampleI++] = procI;
            }
183
184
185
186
187
188
189
190
191
        }
    }
}


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

205
    switch (mode)
206
207
208
    {
        case NEARESTCELL:
        {
209
            if (samplePatch_.size() && samplePatch_ != "none")
210
211
212
213
214
215
            {
                FatalErrorIn
                (
                    "mappedPatchBase::findSamples(const pointField&,"
                    " labelList&, labelList&, pointField&) const"
                )   << "No need to supply a patch name when in "
216
                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
217
218
            }

219
            //- Note: face-diagonal decomposition
220
            const indexedOctree<Foam::treeDataCell>& tree = mesh.cellTree();
221
222
223
224
225

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

226
                label cellI = tree.findInside(sample);
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249

                if (cellI == -1)
                {
                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
                    nearest[sampleI].second().second() = Pstream::myProcNo();
                }
                else
                {
                    const point& cc = mesh.cellCentres()[cellI];

                    nearest[sampleI].first() = pointIndexHit
                    (
                        true,
                        cc,
                        cellI
                    );
                    nearest[sampleI].second().first() = magSqr(cc-sample);
                    nearest[sampleI].second().second() = Pstream::myProcNo();
                }
            }
            break;
        }

250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
        case NEARESTONLYCELL:
        {
            if (samplePatch_.size() && samplePatch_ != "none")
            {
                FatalErrorIn
                (
                    "mappedPatchBase::findSamples(const pointField&,"
                    " labelList&, labelList&, pointField&) const"
                )   << "No need to supply a patch name when in "
                    << 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;
        }

280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
        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,
304
                        1e-4
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
344
345
346
347
348
349
350
351
352
353
                    )
                );
                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;
        }

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
415
416
417
418
419
420
421
422
423
424
        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;
        }

425
426
        case NEARESTFACE:
        {
427
            if (samplePatch().size() && samplePatch() != "none")
428
429
430
431
432
433
            {
                FatalErrorIn
                (
                    "mappedPatchBase::findSamples(const pointField&,"
                    " labelList&, labelList&, pointField&) const"
                )   << "No need to supply a patch name when in "
434
                    << sampleModeNames_[mode] << " mode." << exit(FatalError);
435
436
            }

437
438
439
            //- Note: face-diagonal decomposition
            const meshSearchMeshObject& meshSearchEngine =
                meshSearchMeshObject::New(mesh);
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482

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

                label faceI = meshSearchEngine.findNearestFace(sample);

                if (faceI == -1)
                {
                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
                    nearest[sampleI].second().second() = Pstream::myProcNo();
                }
                else
                {
                    const point& fc = mesh.faceCentres()[faceI];

                    nearest[sampleI].first() = pointIndexHit
                    (
                        true,
                        fc,
                        faceI
                    );
                    nearest[sampleI].second().first() = magSqr(fc-sample);
                    nearest[sampleI].second().second() = Pstream::myProcNo();
                }
            }
            break;
        }

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

        default:
        {
            FatalErrorIn("mappedPatchBase::findSamples(..)")
                << "problem." << abort(FatalError);
        }
    }


483
    // Find nearest. Combine on master.
484
485
486
    Pstream::listCombineGather(nearest, nearestEqOp());
    Pstream::listCombineScatter(nearest);

487

488
489
    if (debug)
    {
490
        Info<< "mappedPatchBase::findSamples on mesh " << sampleRegion()
491
492
493
494
495
496
497
498
            << " : " << endl;
        forAll(nearest, sampleI)
        {
            label procI = nearest[sampleI].second().second();
            label localI = nearest[sampleI].first().index();

            Info<< "    " << sampleI << " coord:"<< samples[sampleI]
                << " found on processor:" << procI
499
500
501
                << " in local cell/face/point:" << localI
                << " with location:" << nearest[sampleI].first().rawPoint()
                << endl;
502
503
504
505
506
507
508
509
510
511
        }
    }

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

    forAll(nearest, sampleI)
    {
512
513
514
515
516
517
518
519
520
521
522
523
        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();
        }
524
525
526
527
528
529
    }
}


void Foam::mappedPatchBase::calcMapping() const
{
530
    static bool hasWarned = false;
531
532
533
534
535
536
    if (mapPtr_.valid())
    {
        FatalErrorIn("mappedPatchBase::calcMapping() const")
            << "Mapping already calculated" << exit(FatalError);
    }

537
538
539
540
541
    // Get points on face (since cannot use face-centres - might be off
    // face-diagonal decomposed tets.
    tmp<pointField> patchPoints(facePoints(patch_));

    // Get offsetted points
542
    const pointField offsettedPoints(samplePoints(patchPoints()));
543

544
545
    // Do a sanity check - am I sampling my own patch?
    // This only makes sense for a non-zero offset.
546
547
548
    bool sampleMyself =
    (
        mode_ == NEARESTPATCHFACE
549
550
     && sampleRegion() == patch_.boundaryMesh().mesh().name()
     && samplePatch() == patch_.name()
551
552
553
    );

    // Check offset
554
555
556
557
    vectorField d(offsettedPoints-patchPoints());
    bool coincident = (gAverage(mag(d)) <= ROOTVSMALL);

    if (sampleMyself && coincident)
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
    {
        WarningIn
        (
            "mappedPatchBase::mappedPatchBase\n"
            "(\n"
            "    const polyPatch& pp,\n"
            "    const word& sampleRegion,\n"
            "    const sampleMode mode,\n"
            "    const word& samplePatch,\n"
            "    const vector& offset\n"
            ")\n"
        )   << "Invalid offset " << d << endl
            << "Offset is the vector added to the patch face centres to"
            << " find the patch face supplying the data." << endl
            << "Setting it to " << d
            << " on the same patch, on the same region"
            << " will find the faces themselves which does not make sense"
            << " for anything but testing." << endl
            << "patch_:" << patch_.name() << endl
577
            << "sampleRegion_:" << sampleRegion() << endl
578
            << "mode_:" << sampleModeNames_[mode_] << endl
579
            << "samplePatch_:" << samplePatch() << endl
580
581
582
583
584
585
586
587
            << "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;
588
589
590
591
592
593
594
595
596
    collectSamples
    (
        patchPoints,
        samples,
        patchFaceProcs,
        patchFaces,
        patchFc
    );

597
598
599
600
    // Find processor and cell/face samples are in and actual location.
    labelList sampleProcs;
    labelList sampleIndices;
    pointField sampleLocations;
601
    findSamples(mode_, samples, sampleProcs, sampleIndices, sampleLocations);
602

603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
    // 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)
            {
                WarningIn
                (
                    "mappedPatchBase::mappedPatchBase\n"
                    "(\n"
                    "    const polyPatch& pp,\n"
                    "    const word& sampleRegion,\n"
                    "    const sampleMode mode,\n"
                    "    const word& samplePatch,\n"
                    "    const vector& offset\n"
                    ")\n"
                )   << "Did not find " << nNotFound
                    << " out of " << sampleProcs.size() << " total samples."
                    << " Sampling these on owner cell centre instead." << endl
634
                    << "On patch " << patch_.name()
635
                    << " on region " << sampleRegion()
636
                    << " in mode " << sampleModeNames_[mode_] << endl
637
638
                    << "with offset mode " << offsetModeNames_[offsetMode_]
                    << ". Suppressing further warnings from " << type() << endl;
639
640
641
642

                hasWarned = true;
            }

643
644
645
            // Collect the samples that cannot be found
            DynamicList<label> subMap;
            DynamicField<point> subSamples;
646
647
648
649
650

            forAll(sampleProcs, sampleI)
            {
                if (sampleProcs[sampleI] == -1)
                {
651
652
                    subMap.append(sampleI);
                    subSamples.append(samples[sampleI]);
653
654
655
                }
            }

656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
            // 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;
673
674
675
        }
    }

676
677
678
679
680
681
682
683
684
685
686
687
688
    // 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
689
    //        << "Found data in region " << sampleRegion()
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
    //        << " 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"
707
708
            << " sampled cell/faceCentres/points to file " << str.name()
            << endl;
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
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

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

    forAll(subMap, procI)
    {
        subMap[procI] = UIndirectList<label>
        (
            sampleIndices,
            subMap[procI]
        );
        constructMap[procI] = UIndirectList<label>
        (
            patchFaces,
            constructMap[procI]
        );

        //if (debug)
        //{
        //    Pout<< "To proc:" << procI << " sending values of cells/faces:"
        //        << subMap[procI] << endl;
        //    Pout<< "From proc:" << procI
        //        << " receiving values of patch faces:"
        //        << constructMap[procI] << endl;
        //}
    }

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

    if (debug)
    {
        // Check that all elements get a value.
        PackedBoolList used(patch_.size());
        forAll(constructMap, procI)
        {
            const labelList& map = constructMap[procI];

            forAll(map, i)
            {
                label faceI = map[i];

                if (used[faceI] == 0)
                {
                    used[faceI] = 1;
                }
                else
                {
                    FatalErrorIn("mappedPatchBase::calcMapping() const")
                        << "On patch " << patch_.name()
                        << " patchface " << faceI
                        << " is assigned to more than once."
                        << abort(FatalError);
                }
            }
        }
        forAll(used, faceI)
        {
            if (used[faceI] == 0)
            {
                FatalErrorIn("mappedPatchBase::calcMapping() const")
                    << "On patch " << patch_.name()
                    << " patchface " << faceI
                    << " is never assigned to."
                    << abort(FatalError);
            }
        }
    }
}


const Foam::autoPtr<Foam::searchableSurface>& Foam::mappedPatchBase::surfPtr()
const
{
801
    const word surfType(surfDict_.lookupOrDefault<word>("type", "none"));
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839

    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())
    {
        FatalErrorIn("mappedPatchBase::calcAMI() const")
            << "AMI already calculated" << exit(FatalError);
    }

    AMIPtr_.clear();

840

841
    const polyPatch& nbr = samplePolyPatch();
842

843
844
    // Transform neighbour patch to local system
    pointField nbrPoints(samplePoints(nbr.localPoints()));
845

846
847
848
    primitivePatch nbrPatch0
    (
        SubList<face>
849
        (
850
851
852
853
854
855
856
857
858
859
860
            nbr.localFaces(),
            nbr.size()
        ),
        nbrPoints
    );


    if (debug)
    {
        OFstream os(patch_.name() + "_neighbourPatch-org.obj");
        meshTools::writeOBJ(os, samplePolyPatch().localFaces(), nbrPoints);
861
862
863
864
865
866
867

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

869
870
871
872
873
874
    // Construct/apply AMI interpolation to determine addressing and weights
    AMIPtr_.reset
    (
        new AMIPatchToPatchInterpolation
        (
            patch_,
875
            nbrPatch0,
876
            surfPtr(),
877
            faceAreaIntersect::tmMesh,
878
            true,
andy's avatar
andy committed
879
            AMIPatchToPatchInterpolation::imFaceAreaWeight,
880
            -1,
881
            AMIReverse_
882
883
884
885
886
        )
    );
}


887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
// 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());
    pointField& fld = tfld();

    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)
                {
                    FatalIOErrorIn
                    (
                        "mappedPatchBase::readListOrField"
                        "(const word& keyword, const dictionary&, const label)",
                        dict
                    )   << "size " << fld.size()
                        << " is not equal to the given value of " << size
                        << exit(FatalIOError);
                }
            }
            else
            {
                FatalIOErrorIn
                (
                    "mappedPatchBase::readListOrField"
                    "(const word& keyword, const dictionary&, const label)",
                    dict
                )   << "expected keyword 'uniform' or 'nonuniform', found "
                    << firstToken.wordToken()
                    << exit(FatalIOError);
            }
        }
        else
        {
            if (is.version() == 2.0)
            {
                IOWarningIn
                (
                    "mappedPatchBase::readListOrField"
                    "(const word& keyword, const dictionary&, const label)",
                    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;
}


963
964
965
966
967
968
969
970
971
972
// * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //

Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp
)
:
    patch_(pp),
    sampleRegion_(patch_.boundaryMesh().mesh().name()),
    mode_(NEARESTPATCHFACE),
973
974
    samplePatch_(""),
    coupleGroup_(),
975
976
977
978
979
    offsetMode_(UNIFORM),
    offset_(vector::zero),
    offsets_(pp.size(), offset_),
    distance_(0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
980
981
    mapPtr_(NULL),
    AMIPtr_(NULL),
982
    AMIReverse_(false),
983
    surfPtr_(NULL),
984
    surfDict_(fileName("surface"))
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
{}


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),
1001
    coupleGroup_(),
1002
1003
1004
1005
1006
    offsetMode_(NONUNIFORM),
    offset_(vector::zero),
    offsets_(offsets),
    distance_(0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
1007
1008
    mapPtr_(NULL),
    AMIPtr_(NULL),
1009
    AMIReverse_(false),
1010
    surfPtr_(NULL),
1011
    surfDict_(fileName("surface"))
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
{}


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),
1028
    coupleGroup_(),
1029
1030
1031
1032
1033
    offsetMode_(UNIFORM),
    offset_(offset),
    offsets_(0),
    distance_(0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
1034
1035
    mapPtr_(NULL),
    AMIPtr_(NULL),
1036
    AMIReverse_(false),
1037
    surfPtr_(NULL),
1038
    surfDict_(fileName("surface"))
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
{}


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),
1055
    coupleGroup_(),
1056
1057
1058
1059
1060
    offsetMode_(NORMAL),
    offset_(vector::zero),
    offsets_(0),
    distance_(distance),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
1061
1062
    mapPtr_(NULL),
    AMIPtr_(NULL),
1063
    AMIReverse_(false),
1064
    surfPtr_(NULL),
1065
    surfDict_(fileName("surface"))
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
    const dictionary& dict
)
:
    patch_(pp),
1076
    sampleRegion_(dict.lookupOrDefault<word>("sampleRegion", "")),
1077
    mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
1078
1079
    samplePatch_(dict.lookupOrDefault<word>("samplePatch", "")),
    coupleGroup_(dict),
1080
1081
1082
1083
1084
1085
1086
    offsetMode_(UNIFORM),
    offset_(vector::zero),
    offsets_(0),
    distance_(0.0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
    mapPtr_(NULL),
    AMIPtr_(NULL),
1087
    AMIReverse_(dict.lookupOrDefault<bool>("flipNormals", false)),
1088
1089
1090
    surfPtr_(NULL),
    surfDict_(dict.subOrEmptyDict("surface"))
{
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
    if (!coupleGroup_.valid())
    {
        if (sampleRegion_.empty())
        {
            // If no coupleGroup and no sampleRegion assume local region
            sampleRegion_ = patch_.boundaryMesh().mesh().name();
            sameRegion_ = true;
        }
    }

1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
    if (dict.found("offsetMode"))
    {
        offsetMode_ = offsetModeNames_.read(dict.lookup("offsetMode"));

        switch(offsetMode_)
        {
            case UNIFORM:
            {
                offset_ = point(dict.lookup("offset"));
            }
            break;

            case NONUNIFORM:
            {
1115
1116
                //offsets_ = pointField(dict.lookup("offsets"));
                offsets_ = readListOrField("offsets", dict, patch_.size());
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
            }
            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;
1135
1136
        //offsets_ = pointField(dict.lookup("offsets"));
        offsets_ = readListOrField("offsets", dict, patch_.size());
1137
    }
1138
    else if (mode_ != NEARESTPATCHFACE && mode_ != NEARESTPATCHFACEAMI)
1139
1140
1141
1142
1143
    {
        FatalIOErrorIn
        (
            "mappedPatchBase::mappedPatchBase\n"
            "(\n"
1144
1145
            "    const polyPatch&,\n"
            "    const dictionary&\n"
1146
1147
1148
1149
1150
1151
1152
1153
1154
            ")\n",
            dict
        )   << "Please supply the offsetMode as one of "
            << NamedEnum<offsetMode, 3>::words()
            << exit(FatalIOError);
    }
}


1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
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),
    offset_(vector::zero),
    offsets_(0),
    distance_(0.0),
    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
    mapPtr_(NULL),
    AMIPtr_(NULL),
    AMIReverse_(dict.lookupOrDefault<bool>("flipNormals", false)),
    surfPtr_(NULL),
    surfDict_(dict.subOrEmptyDict("surface"))
{
    if (mode != NEARESTPATCHFACE && mode != NEARESTPATCHFACEAMI)
    {
        FatalIOErrorIn
        (
            "mappedPatchBase::mappedPatchBase\n"
            "(\n"
            "    const polyPatch&,\n"
            "    const sampleMode,\n"
            "    const dictionary&\n"
            ")\n",
            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;
        }
    }
}


1209
1210
1211
Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
1212
    const mappedPatchBase& mpb
1213
1214
1215
)
:
    patch_(pp),
1216
1217
1218
    sampleRegion_(mpb.sampleRegion_),
    mode_(mpb.mode_),
    samplePatch_(mpb.samplePatch_),
1219
    coupleGroup_(mpb.coupleGroup_),
1220
1221
1222
1223
1224
    offsetMode_(mpb.offsetMode_),
    offset_(mpb.offset_),
    offsets_(mpb.offsets_),
    distance_(mpb.distance_),
    sameRegion_(mpb.sameRegion_),
1225
1226
    mapPtr_(NULL),
    AMIPtr_(NULL),
1227
    AMIReverse_(mpb.AMIReverse_),
1228
    surfPtr_(NULL),
1229
    surfDict_(mpb.surfDict_)
1230
1231
1232
1233
1234
1235
{}


Foam::mappedPatchBase::mappedPatchBase
(
    const polyPatch& pp,
1236
    const mappedPatchBase& mpb,
1237
1238
1239
1240
    const labelUList& mapAddressing
)
:
    patch_(pp),
1241
1242
1243
    sampleRegion_(mpb.sampleRegion_),
    mode_(mpb.mode_),
    samplePatch_(mpb.samplePatch_),
1244
    coupleGroup_(mpb.coupleGroup_),
1245
1246
    offsetMode_(mpb.offsetMode_),
    offset_(mpb.offset_),
1247
1248
1249
1250
1251
1252
    offsets_
    (
        offsetMode_ == NONUNIFORM
      ? vectorField(mpb.offsets_, mapAddressing)
      : vectorField(0)
    ),
1253
1254
    distance_(mpb.distance_),
    sameRegion_(mpb.sameRegion_),
1255
1256
    mapPtr_(NULL),
    AMIPtr_(NULL),
1257
    AMIReverse_(mpb.AMIReverse_),
1258
    surfPtr_(NULL),
1259
    surfDict_(mpb.surfDict_)
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
{}


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

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


void Foam::mappedPatchBase::clearOut()
{
    mapPtr_.clear();
1274
1275
    AMIPtr_.clear();
    surfPtr_.clear();
1276
1277
1278
1279
1280
1281
1282
1283
1284
}


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

const Foam::polyMesh& Foam::mappedPatchBase::sampleMesh() const
{
    return patch_.boundaryMesh().mesh().time().lookupObject<polyMesh>
    (
1285
        sampleRegion()
1286
1287
1288
1289
1290
1291
1292
1293
    );
}


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

1294
    const label patchI = nbrMesh.boundaryMesh().findPatchID(samplePatch());
1295
1296
1297

    if (patchI == -1)
    {
1298
        FatalErrorIn("mappedPatchBase::samplePolyPatch()")
1299
            << "Cannot find patch " << samplePatch()
1300
1301
1302
1303
1304
1305
1306
1307
1308
            << " in region " << sampleRegion_ << endl
            << "Valid patches are " << nbrMesh.boundaryMesh().names()
            << exit(FatalError);
    }

    return nbrMesh.boundaryMesh()[patchI];
}


1309
1310
1311
1312
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints
(
    const pointField& fc
) const
1313
{
1314
    tmp<pointField> tfld(new pointField(fc));
1315
1316
    pointField& fld = tfld();

1317
    switch (offsetMode_)
1318
1319
1320
1321
    {
        case UNIFORM:
        {
            fld += offset_;
1322
            break;
1323
1324
1325
1326
        }
        case NONUNIFORM:
        {
            fld += offsets_;
1327
            break;
1328
1329
1330
1331
1332
1333
1334
1335
        }
        case NORMAL:
        {
            // Get outwards pointing normal
            vectorField n(patch_.faceAreas());
            n /= mag(n);

            fld += distance_*n;
1336
            break;
1337
1338
1339
1340
1341
1342
1343
        }
    }

    return tfld;
}


1344
1345
1346
1347
1348
1349
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints() const
{
    return samplePoints(facePoints(patch_));
}


1350
1351
1352
1353
Foam::pointIndexHit Foam::mappedPatchBase::facePoint
(
    const polyMesh& mesh,
    const label faceI,
1354
    const polyMesh::cellDecomposition decompMode
1355
1356
1357
1358
1359
1360
)
{
    const point& fc = mesh.faceCentres()[faceI];

    switch (decompMode)
    {
1361
1362
        case polyMesh::FACE_PLANES:
        case polyMesh::FACE_CENTRE_TRIS:
1363
1364
1365
1366
1367
1368
1369
        {
            // For both decompositions the face centre is guaranteed to be
            // on the face
            return pointIndexHit(true, fc, faceI);
        }
        break;

1370
1371
        case polyMesh::FACE_DIAG_TRIS:
        case polyMesh::CELL_TETS:
1372
        {
1373
            // Find the intersection of a ray from face centre to cell centre
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
            // Find intersection of (face-centre-decomposition) centre to
            // cell-centre with face-diagonal-decomposition triangles.

            const pointField& p = mesh.points();
            const face& f = mesh.faces()[faceI];

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

            label cellI = mesh.faceOwner()[faceI];
            const point& cc = mesh.cellCentres()[cellI];
            vector d = fc-cc;

            const label fp0 = mesh.tetBasePtIs()[faceI];
            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:
        {
            FatalErrorIn("mappedPatchBase::facePoint()")
                << "problem" << abort(FatalError);
            return pointIndexHit();
        }
    }
}


1431
1432
1433
1434
void Foam::mappedPatchBase::write(Ostream& os) const
{
    os.writeKeyword("sampleMode") << sampleModeNames_[mode_]
        << token::END_STATEMENT << nl;
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
    if (!sampleRegion_.empty())
    {
        os.writeKeyword("sampleRegion") << sampleRegion_
            << token::END_STATEMENT << nl;
    }
    if (!samplePatch_.empty())
    {
        os.writeKeyword("samplePatch") << samplePatch_
            << token::END_STATEMENT << nl;
    }
1445
    coupleGroup_.write(os);
1446

1447
1448
1449
1450
1451
1452
    if
    (
        offsetMode_ == UNIFORM
     && offset_ == vector::zero
     && (mode_ == NEARESTPATCHFACE || mode_ == NEARESTPATCHFACEAMI)
    )
1453
    {
1454
        // Collocated mode. No need to write offset data
1455
    }
1456
    else
1457
    {
1458
1459
1460
1461
        os.writeKeyword("offsetMode") << offsetModeNames_[offsetMode_]
            << token::END_STATEMENT << nl;

        switch (offsetMode_)
1462
        {
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
            case UNIFORM:
            {
                os.writeKeyword("offset") << offset_ << token::END_STATEMENT
                    << nl;
                break;
            }
            case NONUNIFORM:
            {
                offsets_.writeEntry("offsets", os);
                break;
            }
            case NORMAL:
            {
                os.writeKeyword("distance") << distance_ << token::END_STATEMENT
                    << nl;
                break;
            }
1480
1481
        }

1482
        if (mode_ == NEARESTPATCHFACEAMI)
1483
        {
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
            if (AMIReverse_)
            {
                os.writeKeyword("flipNormals") << AMIReverse_
                    << token::END_STATEMENT << nl;
            }

            if (!surfDict_.empty())
            {
                os.writeKeyword(surfDict_.dictName());
                os  << surfDict_;
            }
1495
        }
1496
    }
1497
1498
1499
1500
}


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