mappedPatchBase.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2015-2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "mappedPatchBase.H"
31#include "ListListOps.H"
34#include "meshTools.H"
35#include "OFstream.H"
36#include "Random.H"
37#include "treeDataFace.H"
38#include "treeDataPoint.H"
39#include "indexedOctree.H"
40#include "polyMesh.H"
41#include "polyPatch.H"
42#include "Time.H"
43#include "mapDistribute.H"
44#include "SubField.H"
45#include "triPointRef.H"
46#include "syncTools.H"
47#include "treeDataCell.H"
48#include "DynamicField.H"
49#include "faceAreaWeightAMI.H"
50#include "OTstream.H"
51
52// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
53
54namespace Foam
55{
57}
58
59
60const Foam::Enum
61<
63>
65({
66 { sampleMode::NEARESTCELL, "nearestCell" },
67 { sampleMode::NEARESTPATCHFACE, "nearestPatchFace" },
68 { sampleMode::NEARESTPATCHFACEAMI, "nearestPatchFaceAMI" },
69 { sampleMode::NEARESTPATCHPOINT, "nearestPatchPoint" },
70 { sampleMode::NEARESTFACE, "nearestFace" },
71 { sampleMode::NEARESTONLYCELL, "nearestOnlyCell" },
72});
73
74
75const Foam::Enum
76<
78>
80({
81 { offsetMode::UNIFORM, "uniform" },
82 { offsetMode::NONUNIFORM, "nonuniform" },
83 { offsetMode::NORMAL, "normal" },
84});
85
86
87// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
88
90(
91 const dictionary& dict
92)
93{
94 if (dict.found("sampleDatabase"))
95 {
96 if (dict.get<bool>("sampleDatabase"))
97 {
99 (
101 (
102 "sampleDatabasePath",
104 )
105 );
106 }
107 }
108 else if (dict.found("sampleDatabasePath"))
109 {
110 return autoPtr<fileName>::New(dict.get<fileName>("sampleDatabasePath"));
111 }
112
113 return nullptr;
114}
115
116
118{
119 if (sameWorld())
120 {
121 return true;
122 }
123
124 const Time& runTime = patch_.boundaryMesh().mesh().time();
125 return const_cast<multiWorldConnections&>
126 (
128 ).addConnectionByName(sampleWorld_);
129}
130
131
133{
134 if (sameWorld())
135 {
136 return UPstream::worldComm;
137 }
138
139 const Time& runTime = patch_.boundaryMesh().mesh().time();
140 return
141 multiWorldConnections::New(runTime).getCommByName(sampleWorld_);
142}
143
144
146(
147 const polyPatch& pp
148) const
149{
150 const polyMesh& mesh = pp.boundaryMesh().mesh();
151
152 // Force construction of min-tet decomp
153 (void)mesh.tetBasePtIs();
154
155 // Initialise to face-centre
156 auto tfacePoints = tmp<pointField>::New(patch_.size());
157 auto& facePoints = tfacePoints.ref();
158
159 forAll(pp, facei)
160 {
161 facePoints[facei] = facePoint
162 (
163 mesh,
164 pp.start()+facei,
166 ).rawPoint();
167 }
168
169 return tfacePoints;
170}
171
172
174(
175 const label mySampleWorld, // Wanted world
176 const pointField& facePoints,
177 pointField& samples, // All samples
178 labelList& patchFaceWorlds, // Per sample: wanted world
179 labelList& patchFaceProcs, // Per sample: originating processor
180 labelList& patchFaces, // Per sample: originating patchFace index
181 pointField& patchFc // Per sample: originating centre
182) const
183{
185
186 const label myComm = getCommunicator(); // Get or create
187 const label myRank = Pstream::myProcNo(myComm);
188 const label nProcs = Pstream::nProcs(myComm);
189
190 const label oldWarnComm(Pstream::warnComm);
191 Pstream::warnComm = myComm;
192
193 if (debug & 2)
194 {
195 Perr<< "patch: " << patch_.name()
196 << "[rank=" << myRank << " procs=" << nProcs
197 << " comm=" << myComm << "] collect samples" << endl;
198 }
199
200 // Collect all sample points and the faces they come from.
201 {
202 List<pointField> globalFc(nProcs);
203 globalFc[myRank] = facePoints;
204
205 Pstream::allGatherList(globalFc, Pstream::msgType(), myComm);
206
207 // Rework into straight list
208 patchFc = ListListOps::combine<pointField>
209 (
210 globalFc,
212 );
213 }
214
215 {
216 List<pointField> globalSamples(nProcs);
217 globalSamples[myRank] = samplePoints(facePoints);
218 Pstream::allGatherList(globalSamples, Pstream::msgType(), myComm);
219
220 // Rework into straight list
221 samples = ListListOps::combine<pointField>
222 (
223 globalSamples,
225 );
226 }
227
228 {
229 labelListList globalFaces(nProcs);
230 globalFaces[myRank] = identity(patch_.size());
231 // Distribute to all processors
232 Pstream::allGatherList(globalFaces, Pstream::msgType(), myComm);
233
234 patchFaces = ListListOps::combine<labelList>
235 (
236 globalFaces,
238 );
239 }
240
241 {
242 labelList procToWorldIndex
243 (
244 UPstream::listGatherValues<label>(mySampleWorld, myComm)
245 );
246 labelList nPerProc
247 (
248 UPstream::listGatherValues<label>(patch_.size(), myComm)
249 );
250
251 Pstream::broadcasts(myComm, procToWorldIndex, nPerProc);
252
253 patchFaceWorlds.setSize(patchFaces.size());
254 patchFaceProcs.setSize(patchFaces.size());
255
256 label sampleI = 0;
257 forAll(nPerProc, proci)
258 {
259 for (label i = 0; i < nPerProc[proci]; i++)
260 {
261 patchFaceWorlds[sampleI] = procToWorldIndex[proci];
262 patchFaceProcs[sampleI] = proci;
263 sampleI++;
264 }
265 }
266 }
267
268 Pstream::warnComm = oldWarnComm;
269}
270
271
273(
274 const sampleMode mode,
275
276 const label mySampleWorld, // local world to sample == my own world
277 const word& sampleRegion, // local region to sample
278 const word& samplePatch, // local patch to sample
279
280 const pointField& samples,
281 List<nearInfoWorld>& nearest
282) const
283{
285
286 const label myComm = getCommunicator(); // Get or create
287
288 // Find the local cell containing the samples
289 const label myRank = Pstream::myProcNo(myComm);
290
291 // Lookup the correct region
292 const polyMesh& mesh = lookupMesh(sampleRegion);
293
294 // All the info for nearest. Construct to miss
295 nearest.setSize(samples.size());
296 nearInfoWorld miss;
297 {
298 miss.first().second() = Tuple2<scalar, label>(Foam::sqr(GREAT), -1);
299 miss.second() = -1; // set world to be ignored
300 }
301 nearest = miss;
302
303 switch (mode)
304 {
305 case NEARESTCELL:
306 {
307 if (samplePatch.size() && samplePatch != "none")
308 {
310 << "No need to supply a patch name when in "
311 << sampleModeNames_[mode] << " mode." << exit(FatalError);
312 }
313
314 //- Note: face-diagonal decomposition
316
317 forAll(samples, sampleI)
318 {
319 const point& sample = samples[sampleI];
320 nearInfoWorld& near = nearest[sampleI];
321
322 label celli = tree.findInside(sample);
323
324 if (celli == -1)
325 {
326 near.first().second().first() = Foam::sqr(GREAT);
327 near.first().second().second() = myRank;
328 near.second() = mySampleWorld;
329 }
330 else
331 {
332 const point& cc = mesh.cellCentres()[celli];
333
334 near.first().first() = pointIndexHit
335 (
336 true,
337 cc,
338 celli
339 );
340 near.first().second().first() = magSqr(cc-sample);
341 near.first().second().second() = myRank;
342 near.second() = mySampleWorld;
343 }
344 }
345 break;
346 }
347
348 case NEARESTONLYCELL:
349 {
350 if (samplePatch.size() && samplePatch != "none")
351 {
353 << "No need to supply a patch name when in "
354 << sampleModeNames_[mode] << " mode." << exit(FatalError);
355 }
356
357 //- Note: face-diagonal decomposition
359
360 forAll(samples, sampleI)
361 {
362 const point& sample = samples[sampleI];
363 nearInfoWorld& near = nearest[sampleI];
364
365 near.first().first() = tree.findNearest(sample, sqr(GREAT));
366 near.first().second().first() = magSqr
367 (
368 near.first().first().hitPoint()
369 -sample
370 );
371 near.first().second().second() = myRank;
372 near.second() = mySampleWorld;
373 }
374 break;
375 }
376
377 case NEARESTPATCHFACE:
378 {
379 Random rndGen(123456);
380
381 const polyPatch& pp = lookupPatch(sampleRegion, samplePatch);
382
383 if (pp.empty())
384 {
385 forAll(samples, sampleI)
386 {
387 nearInfoWorld& near = nearest[sampleI];
388 near.first().second().first() = Foam::sqr(GREAT);
389 near.first().second().second() = myRank;
390 near.second() = mySampleWorld;
391 }
392 }
393 else
394 {
395 treeBoundBox patchBb
396 (
398 (
399 rndGen,
400 1e-4
401 )
402 );
403 patchBb.min() -= point::uniform(ROOTVSMALL);
404 patchBb.max() += point::uniform(ROOTVSMALL);
405
406 indexedOctree<treeDataFace> boundaryTree
407 (
408 treeDataFace // all information needed to search faces
409 (
410 false, // do not cache bb
411 mesh,
412 identity(pp.range()) // boundary faces only
413 ),
414 patchBb, // overall search domain
415 8, // maxLevel
416 10, // leafsize
417 3.0 // duplicity
418 );
419
420 forAll(samples, sampleI)
421 {
422 const point& sample = samples[sampleI];
423
424 nearInfoWorld& near = nearest[sampleI];
425 pointIndexHit& nearInfo = near.first().first();
426 nearInfo = boundaryTree.findNearest
427 (
428 sample,
429 magSqr(patchBb.span())
430 );
431
432 if (!nearInfo.hit())
433 {
434 near.first().second().first() = Foam::sqr(GREAT);
435 near.first().second().second() = myRank;
436 near.second() = mySampleWorld;
437 }
438 else
439 {
440 point fc(pp[nearInfo.index()].centre(pp.points()));
441 nearInfo.setPoint(fc);
442 near.first().second().first() = magSqr(fc-sample);
443 near.first().second().second() = myRank;
444 near.second() = mySampleWorld;
445 }
446 }
447 }
448 break;
449 }
450
451 case NEARESTPATCHPOINT:
452 {
453 Random rndGen(123456);
454
455 const polyPatch& pp = lookupPatch(sampleRegion, samplePatch);
456
457 if (pp.empty())
458 {
459 forAll(samples, sampleI)
460 {
461 nearInfoWorld& near = nearest[sampleI];
462 near.first().second().first() = Foam::sqr(GREAT);
463 near.first().second().second() = myRank;
464 near.second() = mySampleWorld;
465 }
466 }
467 else
468 {
469 // patch (local) points
470 treeBoundBox patchBb
471 (
473 (
474 rndGen,
475 1e-4
476 )
477 );
478 patchBb.min() -= point::uniform(ROOTVSMALL);
479 patchBb.max() += point::uniform(ROOTVSMALL);
480
482 (
483 treeDataPoint // all information needed to search faces
484 (
485 mesh.points(),
486 pp.meshPoints() // selection of points to search on
487 ),
488 patchBb, // overall search domain
489 8, // maxLevel
490 10, // leafsize
491 3.0 // duplicity
492 );
493
494 forAll(samples, sampleI)
495 {
496 const point& sample = samples[sampleI];
497
498 nearInfoWorld& near = nearest[sampleI];
499 pointIndexHit& nearInfo = near.first().first();
500 nearInfo = boundaryTree.findNearest
501 (
502 sample,
503 magSqr(patchBb.span())
504 );
505
506 if (!nearInfo.hit())
507 {
508 near.first().second().first() = Foam::sqr(GREAT);
509 near.first().second().second() = myRank;
510 near.second() = mySampleWorld;
511 }
512 else
513 {
514 const point& pt = nearInfo.hitPoint();
515
516 near.first().second().first() = magSqr(pt-sample);
517 near.first().second().second() = myRank;
518 near.second() = mySampleWorld;
519 }
520 }
521 }
522 break;
523 }
524
525 case NEARESTFACE:
526 {
527 if (samplePatch.size() && samplePatch != "none")
528 {
530 << "No need to supply a patch name when in "
531 << sampleModeNames_[mode] << " mode." << exit(FatalError);
532 }
533
534 //- Note: face-diagonal decomposition
535 const meshSearchMeshObject& meshSearchEngine =
537
538 forAll(samples, sampleI)
539 {
540 const point& sample = samples[sampleI];
541 nearInfoWorld& near = nearest[sampleI];
542
543 label facei = meshSearchEngine.findNearestFace(sample);
544
545 if (facei == -1)
546 {
547 near.first().second().first() = Foam::sqr(GREAT);
548 near.first().second().second() = myRank;
549 near.second() = mySampleWorld;
550 }
551 else
552 {
553 const point& fc = mesh.faceCentres()[facei];
554
555 near.first().first() = pointIndexHit(true, fc, facei);
556 near.first().second().first() = magSqr(fc-sample);
557 near.first().second().second() = myRank;
558 near.second() = mySampleWorld;
559 }
560 }
561 break;
562 }
563
564 case NEARESTPATCHFACEAMI:
565 {
566 // nothing to do here
567 return;
568 }
569
570 default:
571 {
573 << "problem." << abort(FatalError);
574 }
575 }
576}
577
578
580(
581 const sampleMode mode,
582 const label myWorld,
583 const pointField& samples,
584 const labelList& wantedWorlds,
585 const labelList& origProcs,
586
587 labelList& sampleProcs,
588 labelList& sampleIndices,
589 pointField& sampleLocations
590) const
591{
593
594 // Find the processor/cell containing the samples. Does not account
595 // for samples being found in two processors.
596
597 const label myComm = getCommunicator(); // Get or create
598 const label myRank = Pstream::myProcNo(myComm);
599 const label nProcs = Pstream::nProcs(myComm);
600
601 const label oldWarnComm(Pstream::warnComm);
602 Pstream::warnComm = myComm;
603
604 wordList samplePatches(nProcs);
605 {
606 samplePatches[myRank] = samplePatch_;
607 Pstream::allGatherList(samplePatches, Pstream::msgType(), myComm);
608 }
609 wordList sampleRegions(nProcs);
610 {
611 sampleRegions[myRank] = sampleRegion_;
612 Pstream::allGatherList(sampleRegions, Pstream::msgType(), myComm);
613 }
614
615
616 // Find all the info for nearest
618 forAll(nearest, samplei)
619 {
620 nearest[samplei].first() = nearInfo
621 (
624 );
625 nearest[samplei].second() = wantedWorlds[samplei];
626 }
627
628
629 // Extract samples to search for locally
630 {
631 DynamicList<label> localMap(samples.size());
632 forAll(wantedWorlds, samplei)
633 {
634 if (wantedWorlds[samplei] == myWorld)
635 {
636 localMap.append(samplei);
637 }
638 }
639
640 if (localMap.size())
641 {
642 pointField localSamples(samples, localMap);
643 labelList localOrigProcs(origProcs, localMap);
644
645 //Assume single patch to sample for now
646 const word localOrigPatch(samplePatches[localOrigProcs[0]]);
647 const word localOrigRegion(sampleRegions[localOrigProcs[0]]);
648 List<nearInfoWorld> localNearest(localSamples.size());
649
650 if (debug)
651 {
652 Pout<< "Searching locally for " << localSamples.size()
653 << " samples on region:" << localOrigRegion
654 << " on patch:" << localOrigPatch << endl;
655 }
656 findLocalSamples
657 (
658 mode,
659 myWorld,
660 localOrigRegion,
661 localOrigPatch,
662 localSamples,
663 localNearest
664 );
665 UIndirectList<nearInfoWorld>(nearest, localMap) = localNearest;
666 }
667 }
668
669
670 // Find nearest - globally consistent
672 (
673 nearest,
676 myComm
677 );
678
679 //if (debug)
680 //{
681 // Pout<< "** After combining:" << endl;
682 // forAll(nearest, samplei)
683 // {
684 // Pout<< " sample:" << samples[samplei]
685 // << " origating from proc:" << origProcs[samplei]
686 // << " to be found on world:" << wantedWorlds[samplei] << nl
687 // << " found on world:" << nearest[samplei].second() << nl
688 // << " found on proc:"
689 // << nearest[samplei].first().second().second() << nl
690 // << " found on patchfacei:"
691 // << nearest[samplei].first().first().index() << nl
692 // << " found at location:"
693 // << nearest[samplei].first().first().rawPoint() << nl;
694 // }
695 // Pout<< endl;
696 //}
697
698 // Convert back into proc+local index
699 sampleProcs.setSize(samples.size());
700 sampleIndices.setSize(samples.size());
701 sampleLocations.setSize(samples.size());
702
703 forAll(nearest, sampleI)
704 {
705 const nearInfo& ni = nearest[sampleI].first();
706
707 if (!ni.first().hit())
708 {
709 sampleProcs[sampleI] = -1;
710 sampleIndices[sampleI] = -1;
711 sampleLocations[sampleI] = vector::max;
712 }
713 else
714 {
715 sampleProcs[sampleI] = ni.second().second();
716 sampleIndices[sampleI] = ni.first().index();
717 sampleLocations[sampleI] = ni.first().hitPoint();
718 }
719 }
720
721 Pstream::warnComm = oldWarnComm;
722}
723
724
726{
727 static bool hasWarned = false;
728 if (mapPtr_)
729 {
731 << "Mapping already calculated" << exit(FatalError);
732 }
733
735
736 const label myComm = getCommunicator(); // Get or create
737
739 //if (sampleDatabase() && !sameWorld())
740 //{
741 // const word syncName("syncObjects");
742 // const polyMesh& mesh = patch_.boundaryMesh().mesh();
743 // const Time& runTime = mesh.time();
744 // const functionObjectList& fos = runTime.functionObjects();
745 //
746 // forAll(fos, i)
747 // {
748 // Pout<< "** FO:" << fos[i].name() << " tpye:" << fos[i].type()
749 // << endl;
750 // }
751 //
752 // if (fos.findObjectID(syncName) == -1)
753 // {
754 // //FatalErrorInFunction
755 // WarningInFunction
756 // << "Did not detect functionObject " << syncName
757 // << ". This is used to synchronise data inbetween worlds"
758 // << endl;
759 // }
760 //}
761
762
763 // Get points on face (since cannot use face-centres - might be off
764 // face-diagonal decomposed tets.
765 tmp<pointField> patchPoints(facePoints(patch_));
766
767 // Get offsetted points
768 const pointField offsettedPoints(samplePoints(patchPoints()));
769
770 // Do a sanity check - am I sampling my own patch?
771 // This only makes sense for a non-zero offset.
772 bool sampleMyself =
773 (
774 mode_ == NEARESTPATCHFACE
775 && sampleWorld() == UPstream::myWorld()
776 && sampleRegion() == patch_.boundaryMesh().mesh().name()
777 && samplePatch() == patch_.name()
778 );
779
780 if (sampleMyself)
781 {
782 // Check offset
783 vectorField d(offsettedPoints-patchPoints());
784 bool coincident = (gAverage(mag(d)) <= ROOTVSMALL);
785
786 if (sampleMyself && coincident)
787 {
789 << "Invalid offset " << d << endl
790 << "Offset is the vector added to the patch face centres to"
791 << " find the patch face supplying the data." << endl
792 << "Setting it to " << d
793 << " on the same patch, on the same region"
794 << " will find the faces themselves which does not make sense"
795 << " for anything but testing." << endl
796 << "patch:" << patch_.name() << endl
797 << "sampleRegion:" << sampleRegion() << endl
798 << "mode:" << sampleModeNames_[mode_] << endl
799 << "samplePatch:" << samplePatch() << endl
800 << "offsetMode:" << offsetModeNames_[offsetMode_] << endl;
801 }
802 }
803
804 // Get local world and world-to-sample in index form
805 const label myWorld = Pstream::myWorldID();
806 const label mySampleWorld = Pstream::allWorlds().find(sampleWorld_);
807
808
809 // Get global list of all samples and the processor and face they come from.
810 pointField samples; // coordinates
811 labelList patchFaceWorlds; // world to sample
812 labelList patchFaceProcs; // originating processor
813 labelList patchFaces; // originating face on processor patch
814 pointField patchFc; // originating face centre
815 collectSamples
816 (
817 mySampleWorld, // world I want to sample
818 patchPoints,
819
820 samples,
821 patchFaceWorlds,
822 patchFaceProcs,
823 patchFaces,
824 patchFc
825 );
826
827 //if (debug)
828 //{
829 // forAll(samples, samplei)
830 // {
831 // Pout<< " sample:" << samples[samplei]
832 // << " origating from proc:" << patchFaceProcs[samplei]
833 // << " face:" << patchFaces[samplei]
834 // << " to be found on world:" << patchFaceWorlds[samplei] << nl;
835 // }
836 //}
837
838 // Find processor and cell/face samples are in and actual location.
839 labelList sampleProcs;
840 labelList sampleIndices;
841 pointField sampleLocations;
842 findSamples
843 (
844 mode_,
845 myWorld, // my world (in index form)
846 samples,
847 patchFaceWorlds,
848 patchFaceProcs,
849
850 sampleProcs,
851 sampleIndices,
852 sampleLocations
853 );
854
855 // Check for samples that were not found. This will only happen for
856 // NEARESTCELL since finds cell containing a location
857 if (mode_ == NEARESTCELL)
858 {
859 label nNotFound = 0;
860 forAll(sampleProcs, sampleI)
861 {
862 if (sampleProcs[sampleI] == -1)
863 {
864 nNotFound++;
865 }
866 }
867 reduce(nNotFound, sumOp<label>(), Pstream::msgType(), myComm);
868
869 if (nNotFound > 0)
870 {
871 if (!hasWarned)
872 {
874 << "Did not find " << nNotFound
875 << " out of " << sampleProcs.size() << " total samples."
876 << " Sampling these on owner cell centre instead." << endl
877 << "On patch " << patch_.name()
878 << " on region " << sampleRegion()
879 << " in mode " << sampleModeNames_[mode_] << endl
880 << "with offset mode " << offsetModeNames_[offsetMode_]
881 << ". Suppressing further warnings from " << type() << endl;
882
883 hasWarned = true;
884 }
885
886 // Collect the samples that cannot be found
887 DynamicList<label> subMap;
888 forAll(sampleProcs, sampleI)
889 {
890 if (sampleProcs[sampleI] == -1)
891 {
892 subMap.append(sampleI);
893 }
894 }
895
896 // And re-search for pure nearest (should not fail)
897 labelList subSampleProcs;
898 labelList subSampleIndices;
899 pointField subSampleLocations;
900 findSamples
901 (
902 NEARESTONLYCELL,
903 myWorld, // my world (in index form)
904
905 pointField(samples, subMap),
906 UIndirectList<label>(patchFaceWorlds, subMap)(),
907 UIndirectList<label>(patchFaceProcs, subMap)(),
908
909 subSampleProcs,
910 subSampleIndices,
911 subSampleLocations
912 );
913
914 // Insert
915 labelUIndList(sampleProcs, subMap) = subSampleProcs;
916 labelUIndList(sampleIndices, subMap) = subSampleIndices;
917 UIndirectList<point>(sampleLocations, subMap) = subSampleLocations;
918 }
919 }
920
921 // Now we have all the data we need:
922 // - where sample originates from (so destination when mapping):
923 // patchFaces, patchFaceProcs.
924 // - cell/face sample is in (so source when mapping)
925 // sampleIndices, sampleProcs.
926
927 if (Pstream::master(myComm))
928 {
929 forAll(samples, i)
930 {
931 if (sampleProcs[i] == -1)
932 {
933 FatalErrorInFunction << "did not find sample "
934 << patchFc[i] << " on patch " << patch_.name()
935 << " on region "
936 << patch_.boundaryMesh().mesh().name()
937 << " on processor " << patchFaceProcs[i]
938 << exit(FatalError);
939 }
940 }
941 }
942
943
944 if (debug && Pstream::master(myComm))
945 {
946 //forAll(samples, i)
947 //{
948 // Pout<< i << " need data in region "
949 // << patch_.boundaryMesh().mesh().name()
950 // << " for proc:" << patchFaceProcs[i]
951 // << " face:" << patchFaces[i]
952 // << " at:" << patchFc[i] << endl
953 // << "Found data in region " << sampleRegion()
954 // << " at proc:" << sampleProcs[i]
955 // << " face:" << sampleIndices[i]
956 // << " at:" << sampleLocations[i]
957 // << nl << endl;
958 //}
959
960 OFstream str
961 (
962 patch_.boundaryMesh().mesh().time().path()
963 / patch_.name()
964 + "_mapped.obj"
965 );
966 Pout<< "Dumping mapping as lines from patch faceCentres to"
967 << " sampled cell/faceCentres/points to file " << str.name()
968 << endl;
969
970 label vertI = 0;
971
972 forAll(patchFc, i)
973 {
974 meshTools::writeOBJ(str, patchFc[i]);
975 vertI++;
976 meshTools::writeOBJ(str, sampleLocations[i]);
977 vertI++;
978 str << "l " << vertI-1 << ' ' << vertI << nl;
979 }
980 }
981
982 // Determine schedule.
983 mapPtr_.reset(new mapDistribute(sampleProcs, patchFaceProcs, myComm));
984
985 // Rework the schedule from indices into samples to cell data to send,
986 // face data to receive.
987
988 labelListList& subMap = mapPtr_().subMap();
989 labelListList& constructMap = mapPtr_().constructMap();
990
991 forAll(subMap, proci)
992 {
993 subMap[proci] = labelUIndList(sampleIndices, subMap[proci]);
994 constructMap[proci] = labelUIndList(patchFaces, constructMap[proci]);
995
996 if (debug)
997 {
998 Pout<< "To proc:" << proci << " sending values of cells/faces:"
999 << subMap[proci] << endl;
1000 Pout<< "From proc:" << proci
1001 << " receiving values of patch faces:"
1002 << constructMap[proci] << endl;
1003 }
1004 }
1005
1006
1007 // Redo constructSize
1008 mapPtr_().constructSize() = patch_.size();
1009
1010 if (debug)
1011 {
1012 // Check that all elements get a value.
1013 bitSet used(patch_.size());
1014 forAll(constructMap, proci)
1015 {
1016 const labelList& map = constructMap[proci];
1017
1018 forAll(map, i)
1019 {
1020 label facei = map[i];
1021
1022 if (used.test(facei))
1023 {
1025 << "On patch " << patch_.name()
1026 << " patchface " << facei
1027 << " is assigned to more than once."
1028 << abort(FatalError);
1029 }
1030 else
1031 {
1032 used.set(facei);
1033 }
1034 }
1035 }
1036 forAll(used, facei)
1037 {
1038 if (!used.test(facei))
1039 {
1041 << "On patch " << patch_.name()
1042 << " patchface " << facei
1043 << " is never assigned to."
1044 << abort(FatalError);
1045 }
1046 }
1047 }
1048}
1049
1050
1052const
1053{
1054 const word surfType(surfDict_.getOrDefault<word>("type", "none"));
1055
1056 if (!surfPtr_ && surfType != "none")
1057 {
1058 word surfName(surfDict_.getOrDefault("name", patch_.name()));
1059
1060 const polyMesh& mesh = patch_.boundaryMesh().mesh();
1061
1062 surfPtr_ =
1064 (
1065 surfType,
1066 IOobject
1067 (
1068 surfName,
1069 mesh.time().constant(),
1070 "triSurface",
1071 mesh,
1074 ),
1075 surfDict_
1076 );
1077 }
1078
1079 return surfPtr_;
1080}
1081
1082
1084{
1085 if (AMIPtr_->upToDate())
1086 {
1088 << "AMI already up-to-date"
1089 << endl;
1090
1091 return;
1092 }
1093
1095
1096 const label myComm = getCommunicator(); // Get or create
1097
1098 const label oldWorldComm(Pstream::worldComm);
1099 const label oldWarnComm(Pstream::warnComm);
1100
1101 // Check if running locally
1102 if (sampleWorld_.empty() || sameWorld())
1103 {
1104 const polyPatch& nbr = samplePolyPatch();
1105
1106 // Transform neighbour patch to local system
1107 const pointField nbrPoints(samplePoints(nbr.localPoints()));
1108
1109 const primitivePatch nbrPatch0
1110 (
1112 (
1113 nbr.localFaces(),
1114 nbr.size()
1115 ),
1116 nbrPoints
1117 );
1118
1119
1120 if (debug)
1121 {
1122 OFstream os(patch_.name() + "_neighbourPatch-org.obj");
1123 meshTools::writeOBJ(os, samplePolyPatch().localFaces(), nbrPoints);
1124
1125 OFstream osN(patch_.name() + "_neighbourPatch-trans.obj");
1126 meshTools::writeOBJ(osN, nbrPatch0, nbrPoints);
1127
1128 OFstream osO(patch_.name() + "_ownerPatch.obj");
1129 meshTools::writeOBJ(osO, patch_.localFaces(), patch_.localPoints());
1130 }
1131
1132 // Construct/apply AMI interpolation to determine addressing and
1133 // weights.
1134
1135 // Change to use inter-world communicator
1136 Pstream::worldComm = myComm;
1138
1139 AMIPtr_->calculate(patch_, nbrPatch0, surfPtr());
1140
1141 Pstream::warnComm = oldWarnComm;
1142 Pstream::worldComm = oldWorldComm;
1143 }
1144 else
1145 {
1146 faceList dummyFaces;
1147 pointField dummyPoints;
1148 const primitivePatch dummyPatch
1149 (
1150 SubList<face>(dummyFaces),
1151 dummyPoints
1152 );
1153
1154 // Change to use inter-world communicator
1155 Pstream::worldComm = myComm;
1157
1158 if (masterWorld())
1159 {
1160 // Construct/apply AMI interpolation to determine addressing
1161 // and weights. Have patch_ for src faces, 0 faces for the
1162 // target side
1163 AMIPtr_->calculate(patch_, dummyPatch, surfPtr());
1164 }
1165 else
1166 {
1167 // Construct/apply AMI interpolation to determine addressing
1168 // and weights. Have 0 faces for src side, patch_ for the tgt
1169 // side
1170 AMIPtr_->calculate(dummyPatch, patch_, surfPtr());
1171 }
1172 // Now the AMI addressing/weights will be from src side (on masterWorld
1173 // processors) to tgt side (on other processors)
1174
1175 Pstream::warnComm = oldWarnComm;
1176 Pstream::worldComm = oldWorldComm;
1177 }
1178}
1179
1180
1182(
1183 const objectRegistry& obr,
1184 const wordList& names,
1185 const label index
1186)
1187{
1188 const objectRegistry& sub = obr.subRegistry(names[index], true);
1189 if (index == names.size()-1)
1190 {
1191 return sub;
1192 }
1193 else
1194 {
1195 return subRegistry(sub, names, index+1);
1196 }
1197}
1198
1199
1200// * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * * * //
1201
1203:
1204 patch_(pp),
1205 sampleWorld_(),
1206 sampleRegion_(patch_.boundaryMesh().mesh().name()),
1207 mode_(NEARESTPATCHFACE),
1208 samplePatch_(),
1209 coupleGroup_(),
1210 sampleDatabasePtr_(),
1211 offsetMode_(UNIFORM),
1212 offset_(Zero),
1213 offsets_(pp.size(), offset_),
1214 distance_(0),
1215 communicator_(-1), // Demand-driven (cached value)
1216 sameRegion_(true),
1217 mapPtr_(nullptr),
1218 AMIReverse_(false),
1219 AMIPtr_(new faceAreaWeightAMI(true, AMIReverse_)),
1220 surfPtr_(nullptr),
1221 surfDict_(fileName("surface"))
1222{
1223 // NOTE: same region, no sample-world. Thus no world-world communication
1224}
1225
1226
1228(
1229 const polyPatch& pp,
1230 const word& sampleRegion,
1231 const sampleMode mode,
1232 const word& samplePatch,
1233 const vectorField& offsets
1234)
1235:
1237 (
1238 pp,
1239 sampleRegion,
1240 mode,
1241 samplePatch,
1242 scalar(0)
1243 )
1244{
1246}
1247
1248
1250(
1251 const polyPatch& pp,
1252 const word& sampleRegion,
1253 const sampleMode mode,
1254 const word& samplePatch,
1255 const vector& uniformOffset
1256)
1257:
1259 (
1260 pp,
1261 sampleRegion,
1262 mode,
1263 samplePatch,
1264 scalar(0)
1265 )
1266{
1267 mappedPatchBase::setOffset(uniformOffset);
1268}
1269
1270
1272(
1273 const polyPatch& pp,
1274 const word& sampleRegion,
1275 const sampleMode mode,
1276 const word& samplePatch,
1277 const scalar normalDistance
1278)
1279:
1280 patch_(pp),
1281 sampleWorld_(),
1282 sampleRegion_(sampleRegion),
1283 mode_(mode),
1284 samplePatch_(samplePatch),
1285 coupleGroup_(),
1286 sampleDatabasePtr_(),
1287 offsetMode_(NORMAL),
1288 offset_(Zero),
1289 offsets_(0),
1290 distance_(normalDistance),
1291 communicator_(-1), // Demand-driven (cached value)
1292 sameRegion_
1293 (
1294 sampleWorld_.empty()
1295 && sampleRegion_ == patch_.boundaryMesh().mesh().name()
1296 ),
1297 mapPtr_(nullptr),
1298 AMIReverse_(false),
1299 AMIPtr_(new faceAreaWeightAMI(true, AMIReverse_)),
1300 surfPtr_(nullptr),
1301 surfDict_(fileName("surface"))
1302{
1304}
1305
1306
1308(
1309 const polyPatch& pp,
1310 const dictionary& dict
1311)
1312:
1313 patch_(pp),
1314 sampleWorld_(dict.getOrDefault<word>("sampleWorld", word::null)),
1315 sampleRegion_(dict.getOrDefault<word>("sampleRegion", word::null)),
1316 mode_(sampleModeNames_.get("sampleMode", dict)),
1317 samplePatch_(dict.getOrDefault<word>("samplePatch", word::null)),
1318 coupleGroup_(dict),
1319 sampleDatabasePtr_(readDatabase(dict)),
1320 offsetMode_(UNIFORM),
1321 offset_(Zero),
1322 offsets_(0),
1323 distance_(0),
1324 communicator_(-1), // Demand-driven (cached value)
1325 sameRegion_
1326 (
1327 sampleWorld_.empty()
1328 && sampleRegion_ == patch_.boundaryMesh().mesh().name()
1329 ),
1330 mapPtr_(nullptr),
1331 AMIReverse_
1332 (
1333 dict.getOrDefault
1334 (
1335 "reverseTarget", // AMIInterpolation uses this keyword
1336 dict.getOrDefault
1337 (
1338 "flipNormals",
1339 false
1340 )
1341 )
1342 ),
1343 AMIPtr_
1344 (
1346 (
1347 dict.getOrDefault("AMIMethod", faceAreaWeightAMI::typeName),
1348 dict,
1349 AMIReverse_
1350 )
1351 ),
1352 surfPtr_(nullptr),
1353 surfDict_(dict.subOrEmptyDict("surface"))
1354{
1356
1357 if (!coupleGroup_.valid())
1358 {
1359 if (sampleWorld_.empty() && sampleRegion_.empty())
1360 {
1361 // If no coupleGroup and no sampleRegion assume local region
1363 sameRegion_ = true;
1364 }
1365 }
1366
1367 if (offsetModeNames_.readIfPresent("offsetMode", dict, offsetMode_))
1368 {
1369 switch (offsetMode_)
1370 {
1371 case UNIFORM:
1372 {
1373 dict.readEntry("offset", offset_);
1374 }
1375 break;
1376
1377 case NONUNIFORM:
1378 {
1379 offsets_ = pointField("offsets", dict, patch_.size());
1380 }
1381 break;
1382
1383 case NORMAL:
1384 {
1385 dict.readEntry("distance", distance_);
1386 }
1387 break;
1388 }
1389 }
1390 else if (dict.readIfPresent("offset", offset_))
1391 {
1393 }
1394 else if (dict.found("offsets"))
1395 {
1397 offsets_ = pointField("offsets", dict, patch_.size());
1398 }
1400 {
1402 << "Please supply the offsetMode as one of "
1404 << exit(FatalIOError);
1405 }
1406}
1407
1408
1410(
1411 const polyPatch& pp,
1412 const sampleMode mode,
1413 const dictionary& dict
1414)
1415:
1416 patch_(pp),
1417 sampleWorld_(dict.getOrDefault<word>("sampleWorld", word::null)),
1418 sampleRegion_(dict.getOrDefault<word>("sampleRegion", word::null)),
1419 mode_(mode),
1420 samplePatch_(dict.getOrDefault<word>("samplePatch", word::null)),
1421 coupleGroup_(dict), //dict.getOrDefault<word>("coupleGroup", word::null)),
1422 sampleDatabasePtr_(readDatabase(dict)),
1423 offsetMode_(UNIFORM),
1424 offset_(Zero),
1425 offsets_(0),
1426 distance_(0),
1427 communicator_(-1), // Demand-driven (cached value)
1428 sameRegion_
1429 (
1430 sampleWorld_.empty()
1431 && sampleRegion_ == patch_.boundaryMesh().mesh().name()
1432 ),
1433 mapPtr_(nullptr),
1434 AMIReverse_
1435 (
1436 dict.getOrDefault
1437 (
1438 "reverseTarget", // AMIInterpolation uses this keyword
1439 dict.getOrDefault
1440 (
1441 "flipNormals",
1442 false
1443 )
1444 )
1445 ),
1446 AMIPtr_
1447 (
1449 (
1450 dict.getOrDefault("AMIMethod", faceAreaWeightAMI::typeName),
1451 dict,
1452 AMIReverse_
1453 )
1454 ),
1455 surfPtr_(nullptr),
1456 surfDict_(dict.subOrEmptyDict("surface"))
1457{
1459
1461 {
1463 << "Construct from sampleMode and dictionary only applicable for "
1464 << " collocated patches in modes "
1467 << exit(FatalIOError);
1468 }
1469
1470 if (!coupleGroup_.valid())
1471 {
1472 if (sampleWorld_.empty() && sampleRegion_.empty())
1473 {
1474 // If no coupleGroup and no sampleRegion assume local region
1476 sameRegion_ = true;
1477 }
1478 }
1479}
1480
1481
1483(
1484 const polyPatch& pp,
1485 const mappedPatchBase& mpb
1486)
1487:
1488 patch_(pp),
1489 sampleWorld_(mpb.sampleWorld_),
1490 sampleRegion_(mpb.sampleRegion_),
1491 mode_(mpb.mode_),
1492 samplePatch_(mpb.samplePatch_),
1493 coupleGroup_(mpb.coupleGroup_),
1494 sampleDatabasePtr_
1495 (
1496 mpb.sampleDatabasePtr_
1497 ? new fileName(mpb.sampleDatabasePtr_())
1498 : nullptr
1499 ),
1500 offsetMode_(mpb.offsetMode_),
1501 offset_(mpb.offset_),
1502 offsets_(mpb.offsets_),
1503 distance_(mpb.distance_),
1504 communicator_(mpb.communicator_),
1505 sameRegion_(mpb.sameRegion_),
1506 mapPtr_(nullptr),
1507 AMIReverse_(mpb.AMIReverse_),
1508 AMIPtr_(mpb.AMIPtr_->clone()),
1509 surfPtr_(nullptr),
1510 surfDict_(mpb.surfDict_)
1511{}
1512
1513
1515(
1516 const polyPatch& pp,
1517 const mappedPatchBase& mpb,
1518 const labelUList& mapAddressing
1519)
1520:
1521 patch_(pp),
1522 sampleWorld_(mpb.sampleWorld_),
1523 sampleRegion_(mpb.sampleRegion_),
1524 mode_(mpb.mode_),
1525 samplePatch_(mpb.samplePatch_),
1526 coupleGroup_(mpb.coupleGroup_),
1527 sampleDatabasePtr_
1528 (
1529 mpb.sampleDatabasePtr_
1530 ? new fileName(mpb.sampleDatabasePtr_())
1531 : nullptr
1532 ),
1533 offsetMode_(mpb.offsetMode_),
1534 offset_(mpb.offset_),
1535 offsets_
1536 (
1537 offsetMode_ == NONUNIFORM
1538 ? vectorField(mpb.offsets_, mapAddressing)
1539 : vectorField()
1540 ),
1541 distance_(mpb.distance_),
1542 communicator_(mpb.communicator_),
1543 sameRegion_(mpb.sameRegion_),
1544 mapPtr_(nullptr),
1545 AMIReverse_(mpb.AMIReverse_),
1546 AMIPtr_(mpb.AMIPtr_->clone()),
1547 surfPtr_(nullptr),
1548 surfDict_(mpb.surfDict_)
1549{}
1550
1551
1552// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1553
1555{
1556 clearOut();
1557}
1558
1559
1561{
1562 mapPtr_.reset(nullptr);
1563 surfPtr_.reset(nullptr);
1564 AMIPtr_->upToDate() = false;
1565}
1566
1567
1568// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1569
1570void Foam::mappedPatchBase::setOffset(const scalar normalDist)
1571{
1572 clearOut();
1573 offsetMode_ = offsetMode::NORMAL;
1574 offset_ = Zero;
1575 offsets_.clear();
1576 distance_ = normalDist;
1577}
1578
1579
1581{
1582 clearOut();
1583 offsetMode_ = offsetMode::UNIFORM;
1584 offset_ = uniformOffset;
1585 offsets_.clear();
1586 distance_ = Zero;
1587}
1588
1589
1591{
1592 clearOut();
1593 offsetMode_ = offsetMode::NONUNIFORM;
1594 offset_ = Zero;
1595 offsets_ = offsets;
1596 distance_ = Zero;
1597}
1598
1599
1601(
1602 const word& sampleRegion
1603) const
1604{
1605 const polyMesh& thisMesh = patch_.boundaryMesh().mesh();
1606 return
1607 (
1608 sampleRegion.empty() || sampleRegion == thisMesh.name()
1609 ? thisMesh
1610 : thisMesh.time().lookupObject<polyMesh>(sampleRegion)
1611 );
1612}
1613
1614
1616(
1617 const word& sampleRegion,
1618 const word& samplePatch
1619) const
1620{
1621 const polyMesh& nbrMesh = lookupMesh(sampleRegion);
1622
1623 const label patchi = nbrMesh.boundaryMesh().findPatchID(samplePatch);
1624
1625 if (patchi == -1)
1626 {
1628 << "Cannot find patch " << samplePatch
1629 << " in region " << sampleRegion_ << endl
1630 << exit(FatalError);
1631 }
1632 return nbrMesh.boundaryMesh()[patchi];
1633}
1634
1635
1637{
1638 if (UPstream::myWorld() != sampleWorld_)
1639 {
1641 << "sampleWorld : " << sampleWorld_
1642 << " is not the current world : " << UPstream::myWorld()
1643 << exit(FatalError);
1644 }
1645 return lookupMesh(sampleRegion());
1646}
1647
1648
1650{
1651 const polyMesh& nbrMesh = sampleMesh();
1652
1653 const label patchi = nbrMesh.boundaryMesh().findPatchID(samplePatch());
1654
1655 if (patchi == -1)
1656 {
1658 << "Cannot find patch " << samplePatch()
1659 << " in region " << sampleRegion_ << endl
1660 << "Valid patches are " << nbrMesh.boundaryMesh().names()
1661 << exit(FatalError);
1662 }
1663
1664 return nbrMesh.boundaryMesh()[patchi];
1665}
1666
1667
1669(
1670 const pointField& fc
1671) const
1672{
1673 auto tfld = tmp<pointField>::New(fc);
1674 auto& fld = tfld.ref();
1675
1676 switch (offsetMode_)
1677 {
1678 case UNIFORM:
1679 {
1680 fld += offset_;
1681 break;
1682 }
1683 case NONUNIFORM:
1684 {
1685 fld += offsets_;
1686 break;
1687 }
1688 case NORMAL:
1689 {
1690 // Get outwards pointing normal
1691 vectorField n(patch_.faceAreas());
1692 n /= mag(n);
1693
1694 fld += distance_*n;
1695 break;
1696 }
1697 }
1698
1699 return tfld;
1700}
1701
1702
1704{
1705 return samplePoints(facePoints(patch_));
1706}
1707
1708
1710(
1711 const polyMesh& mesh,
1712 const label facei,
1713 const polyMesh::cellDecomposition decompMode
1714)
1715{
1716 const point& fc = mesh.faceCentres()[facei];
1717
1718 switch (decompMode)
1719 {
1722 {
1723 // For both decompositions the face centre is guaranteed to be
1724 // on the face
1725 return pointIndexHit(true, fc, facei);
1726 }
1727 break;
1728
1731 {
1732 // Find the intersection of a ray from face centre to cell centre
1733 // Find intersection of (face-centre-decomposition) centre to
1734 // cell-centre with face-diagonal-decomposition triangles.
1735
1736 const pointField& p = mesh.points();
1737 const face& f = mesh.faces()[facei];
1738
1739 if (f.size() <= 3)
1740 {
1741 // Return centre of triangle.
1742 return pointIndexHit(true, fc, 0);
1743 }
1744
1745 const label celli = mesh.faceOwner()[facei];
1746 const point& cc = mesh.cellCentres()[celli];
1747 vector d = fc-cc;
1748
1749 const label fp0 = mesh.tetBasePtIs()[facei];
1750 const point& basePoint = p[f[fp0]];
1751
1752 label fp = f.fcIndex(fp0);
1753 for (label i = 2; i < f.size(); i++)
1754 {
1755 const point& thisPoint = p[f[fp]];
1756 label nextFp = f.fcIndex(fp);
1757 const point& nextPoint = p[f[nextFp]];
1758
1759 const triPointRef tri(basePoint, thisPoint, nextPoint);
1760 pointHit hitInfo = tri.intersection
1761 (
1762 cc,
1763 d,
1765 );
1766
1767 if (hitInfo.hit() && hitInfo.distance() > 0)
1768 {
1769 return pointIndexHit(true, hitInfo.hitPoint(), i-2);
1770 }
1771
1772 fp = nextFp;
1773 }
1774
1775 // Fall-back
1776 return pointIndexHit(false, fc, -1);
1777 }
1778 break;
1779
1780 default:
1781 {
1783 << "problem" << abort(FatalError);
1784 return pointIndexHit();
1785 }
1786 }
1787}
1788
1789
1791(
1792 const objectRegistry& obr,
1793 const fileName& path
1794)
1795{
1796 // Lookup (and create if non-existing) a registry using
1797 // '/' separated path. Like 'mkdir -p'
1798
1799 fileName cleanedPath(path);
1800 cleanedPath.clean(); // Remove unneeded ".."
1801 const wordList names(cleanedPath.components());
1802
1803 if (names.empty())
1804 {
1805 return obr;
1806 }
1807 else
1808 {
1809 return subRegistry(obr, names, 0);
1810 }
1811}
1812
1813
1815(
1816 const fileName& root,
1817 const label proci
1818)
1819{
1820 const word processorName("processor"+Foam::name(proci));
1821 return root/"send"/processorName;
1822}
1823
1824
1826{
1827 return sendPath(sampleDatabasePath(), proci);
1828}
1829
1830
1832(
1833 const fileName& root,
1834 const label proci
1835)
1836{
1837 const word processorName("processor"+Foam::name(proci));
1838 return root/"receive"/processorName;
1839}
1840
1841
1843{
1844 return receivePath(sampleDatabasePath(), proci);
1845}
1846
1847
1849(
1850 const objectRegistry& obr,
1852)
1853{
1854 forAllIters(obr, iter)
1855 {
1856 regIOobject* objPtr = iter.val();
1857 const regIOobject& obj = *objPtr;
1858
1859 if (isA<objectRegistry>(obj))
1860 {
1861 dictionary& d = dict.subDictOrAdd(obj.name());
1862 writeDict(dynamic_cast<const objectRegistry&>(obj), d);
1863 }
1864 else if
1865 (
1866 writeIOField<scalar>(obj, dict)
1867 || writeIOField<vector>(obj, dict)
1868 || writeIOField<sphericalTensor>(obj, dict)
1869 || writeIOField<symmTensor>(obj, dict)
1870 || writeIOField<tensor>(obj, dict)
1871 )
1872 {
1873 // IOField specialisation
1874 }
1875 else
1876 {
1877 // Not tested. No way of retrieving data anyway ...
1878 OTstream os;
1879 obj.writeData(os);
1880
1881 primitiveEntry* pePtr = new primitiveEntry(obj.name(), os.tokens());
1882 dict.add(pePtr);
1883 }
1884 }
1885}
1886
1887
1889{
1890 // Reverse of writeDict - reads fields from dictionary into objectRegistry
1891 for (const entry& e : d)
1892 {
1893 if (e.isDict())
1894 {
1895 // Add sub registry
1896 objectRegistry& sub = const_cast<objectRegistry&>
1897 (
1898 obr.subRegistry(e.keyword(), true)
1899 );
1900
1901 readDict(e.dict(), sub);
1902 }
1903 else
1904 {
1905 ITstream& is = e.stream();
1906 token tok(is);
1907
1908 if
1909 (
1910 constructIOField<scalar>(e.keyword(), tok, is, obr)
1911 || constructIOField<vector>(e.keyword(), tok, is, obr)
1912 || constructIOField<sphericalTensor>(e.keyword(), tok, is, obr)
1913 || constructIOField<symmTensor>(e.keyword(), tok, is, obr)
1914 || constructIOField<tensor>(e.keyword(), tok, is, obr)
1915 )
1916 {
1917 // Done storing field
1918 }
1919 else
1920 {
1921 FatalErrorInFunction << "Unsupported type " << e.keyword()
1922 << exit(FatalError);
1923 }
1924 }
1925 }
1926}
1927
1928
1930{
1931 os.writeEntry("sampleMode", sampleModeNames_[mode_]);
1932 os.writeEntryIfDifferent<word>("sampleWorld", word::null, sampleWorld_);
1933 os.writeEntryIfDifferent<word>("sampleRegion", word::null, sampleRegion_);
1934 os.writeEntryIfDifferent<word>("samplePatch", word::null, samplePatch_);
1935
1936 if (sampleDatabasePtr_)
1937 {
1938 os.writeEntry("sampleDatabase", Switch::name(true));
1939 // Write database path if differing
1941 (
1942 "sampleDatabasePath",
1944 sampleDatabasePtr_()
1945 );
1946 }
1947 coupleGroup_.write(os);
1948
1949 if
1950 (
1951 offsetMode_ == UNIFORM
1952 && offset_ == vector::zero
1953 && (mode_ == NEARESTPATCHFACE || mode_ == NEARESTPATCHFACEAMI)
1954 )
1955 {
1956 // Collocated mode. No need to write offset data
1957 }
1958 else
1959 {
1960 os.writeEntry("offsetMode", offsetModeNames_[offsetMode_]);
1961
1962 switch (offsetMode_)
1963 {
1964 case UNIFORM:
1965 {
1966 os.writeEntry("offset", offset_);
1967 break;
1968 }
1969 case NONUNIFORM:
1970 {
1971 offsets_.writeEntry("offsets", os);
1972 break;
1973 }
1974 case NORMAL:
1975 {
1976 os.writeEntry("distance", distance_);
1977 break;
1978 }
1979 }
1980 }
1981
1982 if (mode_ == NEARESTPATCHFACEAMI)
1983 {
1984 if (AMIPtr_)
1985 {
1986 // Use AMI to write itself. Problem: outputs:
1987 // - restartUncoveredSourceFace
1988 // - reverseTarget (instead of flipNormals)
1989 AMIPtr_->write(os);
1990 }
1991 if (!surfDict_.empty())
1992 {
1993 surfDict_.writeEntry(surfDict_.dictName(), os);
1994 }
1995 }
1996}
1997
1998
1999// ************************************************************************* //
label n
Macros for easy insertion into run-time selection tables.
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
Interpolation class dealing with transfer of data between two primitive patches with an arbitrary mes...
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: Enum.H:61
bool readIfPresent(const word &key, const dictionary &dict, EnumType &val) const
Find an entry if present, and assign to T val.
Definition: EnumI.H:132
Minimal example by using system/controlDict.functions:
bool test(const Key &key) const noexcept
Same as found() - return true if key exists in the set.
Definition: HashSet.H:180
bool set(const Key &key)
Same as insert (no value to overwrite)
Definition: HashSet.H:197
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
virtual const fileName & name() const
Return the name of the stream.
Definition: IOstream.C:40
An input stream of tokens.
Definition: ITstream.H:56
void setSize(const label n)
Alias for resize()
Definition: List.H:218
Output to file stream, using an OSstream.
Definition: OFstream.H:57
virtual const fileName & name() const
Read/write access to the name of the stream.
Definition: OSstream.H:107
virtual const fileName & name() const
Get the name of the stream.
Definition: OSstream.H:107
A simple output token stream that can be used to build token lists. Always UNCOMPRESSED.
Definition: OTstream.H:60
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
Definition: Ostream.H:239
Ostream & writeEntryIfDifferent(const word &key, const T &value1, const T &value2)
Write a keyword/value entry only when the two values differ.
Definition: Ostream.H:251
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
Describes the interaction of a face and a point. It carries the info of a successful hit and (if succ...
Definition: PointHit.H:54
scalar distance() const noexcept
Return distance to hit.
Definition: PointHit.H:139
bool hit() const noexcept
Is there a hit.
Definition: PointHit.H:121
const point_type & hitPoint() const
Return the hit point. Fatal if not hit.
Definition: PointHit.H:145
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
A list of faces which address into the list of points.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
static void allGatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
static void listCombineAllGather(const List< commsStruct > &comms, List< T > &values, const CombineOp &cop, const int tag, const label comm)
After completion all processors have the same data.
static void broadcasts(const label comm, Type &arg1, Args &&... args)
Broadcast multiple items to all processes in communicator.
Random number generator.
Definition: Random.H:60
A List obtained as a section of another List.
Definition: SubList.H:70
const word & constant() const
Return constant name.
Definition: TimePathsI.H:96
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
A 2-tuple for storing two objects of dissimilar types. The container is similar in purpose to std::pa...
Definition: Tuple2.H:58
const T1 & first() const noexcept
Return first.
Definition: Tuple2.H:118
const T2 & second() const noexcept
Return second.
Definition: Tuple2.H:130
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
T & first()
Return the first element of the list.
Definition: UListI.H:202
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Definition: UListI.H:427
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label fcIndex(const label i) const noexcept
Definition: UListI.H:60
static label warnComm
Debugging: warn for use of any communicator differing from warnComm.
Definition: UPstream.H:296
static int & msgType() noexcept
Message tag of standard messages.
Definition: UPstream.H:556
static label worldComm
Default communicator (all processors)
Definition: UPstream.H:293
static label myWorldID()
My worldID.
Definition: UPstream.H:495
static const wordList & allWorlds() noexcept
All worlds.
Definition: UPstream.H:483
static const word & myWorld()
My world.
Definition: UPstream.H:501
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:91
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:127
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface....
Definition: boundaryMesh.H:63
bool valid() const noexcept
Is a valid patchGroup (non-empty) name.
const DynamicList< point > & facePoints()
Returns the points of the cutting PLICface.
Definition: cutCellIso.C:173
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
dictionary & subDictOrAdd(const word &keyword, enum keyType::option matchOpt=keyType::REGEX)
Find and return a sub-dictionary for manipulation.
Definition: dictionary.C:500
T lookupOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionary.H:1290
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Search for an entry (const access) with the given keyword.
Definition: dictionaryI.H:87
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
entry * add(entry *entryPtr, bool mergeEntry=false)
Add a new entry.
Definition: dictionary.C:640
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
Particle-size distribution model wherein random samples are drawn from the doubly-truncated uniform p...
Definition: uniform.H:164
A keyword and a list of tokens is an 'entry'.
Definition: entry.H:70
Face area weighted Arbitrary Mesh Interface (AMI) method.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A class for handling file names.
Definition: fileName.H:76
wordList components(const char delim='/') const
Return path components as wordList.
Definition: fileName.C:514
static bool clean(std::string &str)
Definition: fileName.C:199
static const fileName null
An empty fileName.
Definition: fileName.H:102
virtual bool write()
Write the output fields.
Foam::dictionary writeDict() const
Write to dictionary.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
Non-pointer based hierarchical recursive searching.
Definition: indexedOctree.H:74
label findInside(const point &) const
Find shape containing point. Only implemented for certain.
Class containing processor-to-processor mapping information.
Determines a mapping between patch face centres and mesh cell or face centres and processors they're ...
label getWorldCommunicator() const
Get the communicator for the world-world connection.
scalar distance_
Offset distance (normal)
static const objectRegistry & subRegistry(const objectRegistry &obr, const wordList &names, const label index)
const polyPatch & patch_
Patch to sample.
vectorField offsets_
Offset vector (nonuniform)
const vectorField & offsets() const noexcept
Offset vectors (from patch faces to destination mesh objects)
void findSamples(const sampleMode mode, const label myWorldIndex, const pointField &, const labelList &wantedWorlds, const labelList &origProcs, labelList &sampleProcs, labelList &sampleIndices, pointField &sampleLocations) const
Find (global) cells/faces containing samples.
static void readDict(const dictionary &d, objectRegistry &obr)
(recursively) construct and register IOFields from dictionary
sampleMode
Mesh items to sample.
@ NEARESTPATCHFACE
nearest face on selected patch
@ NEARESTPATCHFACEAMI
nearest patch face + AMI interpolation
void findLocalSamples(const sampleMode mode, const label sampleWorld, const word &sampleRegion, const word &samplePatch, const pointField &samplePoints, List< nearInfoWorld > &nearest) const
Find (local) cells/faces containing samples.
void calcAMI() const
Calculate AMI interpolator.
static fileName sendPath(const fileName &root, const label proci)
Helper: return path to store data to be sent to processor i.
static autoPtr< fileName > readDatabase(const dictionary &dict)
Read optional database name from dictionary.
static const Enum< sampleMode > sampleModeNames_
const polyPatch & lookupPatch(const word &sampleRegion, const word &samplePatch) const
Lookup patch.
const sampleMode mode_
What to sample.
const polyMesh & sampleMesh() const
Get the region mesh.
bool sameRegion_
Same region.
word sampleRegion_
Region to sample.
tmp< pointField > samplePoints() const
Get the sample points.
bool addWorldConnection()
Add a world-world connection.
const polyPatch & samplePolyPatch() const
Get the patch on the region.
void collectSamples(const label mySampleWorld, const pointField &facePoints, pointField &samples, labelList &patchFaceWorlds, labelList &patchFaceProcs, labelList &patchFaces, pointField &patchFc) const
const autoPtr< Foam::searchableSurface > & surfPtr() const
Return a pointer to the AMI projection surface.
static fileName receivePath(const fileName &root, const label proci)
offsetMode
How to project face centres.
@ NORMAL
use face normal + distance
@ UNIFORM
single offset vector
@ NONUNIFORM
per-face offset vector
static const Enum< offsetMode > offsetModeNames_
virtual ~mappedPatchBase()
Destructor.
void calcMapping() const
Calculate mapping.
static pointIndexHit facePoint(const polyMesh &, const label facei, const polyMesh::cellDecomposition)
Get a point on the face given a face decomposition method:
sampleMode mode() const noexcept
What to sample.
const coupleGroupIdentifier coupleGroup_
PatchGroup (if in sampleMode NEARESTPATCH*)
void setOffset(const scalar normalDist)
Change to normal offset with given distance.
const polyMesh & lookupMesh(const word &region) const
Lookup mesh.
vector offset_
Offset vector (uniform)
offsetMode offsetMode_
How to obtain samples.
word sampleWorld_
World to sample.
MeshObject wrapper around meshSearch(mesh).
label findNearestFace(const point &location, const label seedFacei=-1, const bool useTreeSearch=true) const
Definition: meshSearch.C:757
Centralized handling of multi-world MPI connections.
Registry of regIOobjects.
const Time & time() const noexcept
Return time registry.
const objectRegistry & subRegistry(const word &name, const bool forceCreate=false, const bool recursive=false) const
Lookup and return a const sub-objectRegistry.
const Type & lookupObject(const word &name, const bool recursive=false) const
static const complex max
complex (VGREAT,VGREAT)
Definition: complex.H:293
label findPatchID(const word &patchName, const bool allowNotFound=true) const
Find patch index given a name, return -1 if not found.
const polyMesh & mesh() const noexcept
Return the mesh reference.
wordList names() const
Return a list of patch names.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
cellDecomposition
Enumeration defining the decomposition of the cell for.
Definition: polyMesh.H:101
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
const labelIOList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:906
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
const indexedOctree< treeDataCell > & cellTree() const
Return the cell search tree.
Definition: polyMesh.C:940
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:364
const polyBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
Definition: polyPatch.C:315
labelRange range() const
Return start/size range of this patch.
Definition: polyPatch.H:370
A keyword and a list of tokens comprise a primitiveEntry. A primitiveEntry can be read,...
const vectorField & faceCentres() const
const vectorField & cellCentres() const
int myProcNo() const noexcept
Return processor number.
regIOobject is an abstract class derived from IOobject to handle automatic object registration with t...
Definition: regIOobject.H:76
virtual bool writeData(Ostream &) const =0
Pure virtual writeData function.
splitCell * master() const
Definition: splitCell.H:113
A class for managing temporary objects.
Definition: tmp.H:65
A token holds an item read from Istream.
Definition: token.H:69
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:89
treeBoundBox extend(Random &rndGen, const scalar s) const
Return slightly wider bounding box.
Encapsulation of data needed to search for faces.
Definition: treeDataFace.H:60
Holds (reference to) pointField. Encapsulation of data needed for octree searches....
Definition: treeDataPoint.H:63
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:80
pointHit intersection(const point &p, const vector &q, const intersection::algorithm alg, const scalar tol=0.0) const
Fast intersection with a ray.
Definition: triangleI.H:439
A class for handling words, derived from Foam::string.
Definition: word.H:68
static const word null
An empty word.
Definition: word.H:80
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
volScalarField & p
dynamicFvMesh & mesh
engineTime & runTime
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:473
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
#define WarningInFunction
Report a warning using Foam::Warning.
#define DebugInFunction
Report an information message using Foam::Info.
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
Namespace for OpenFOAM.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
prefixOSstream Perr
OSstream wrapped stderr (std::cerr) with parallel prefix.
label facePoint(const int facei, const block &block, const label i, const label j)
mode_t mode(const fileName &name, const bool followLink=true)
Return the file mode, normally following symbolic links.
Definition: MSwindows.C:572
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:598
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
errorManip< error > abort(error &err)
Definition: errorManip.H:144
IOerror FatalIOError
Type gAverage(const FieldField< Field, Type > &f)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: IndirectList.H:68
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
dictionary dict
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllIters(container, iter)
Iterate across all elements in the container object.
Definition: stdFoam.H:260
A non-counting (dummy) refCount.
Definition: refCount.H:59
scalarField samples(nIntervals, Zero)
Random rndGen
Definition: createFields.H:23