backgroundMeshDecomposition.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) 2017-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
30#include "meshSearch.H"
33#include "Time.H"
34#include "Random.H"
35#include "pointConversion.H"
36#include "decompositionModel.H"
37
38// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39
40namespace Foam
41{
42 defineTypeNameAndDebug(backgroundMeshDecomposition, 0);
43}
44
45
46// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
47
49(
50 const List<label>& toProc
51)
52{
53 // Determine send map
54 // ~~~~~~~~~~~~~~~~~~
55
56 // 1. Count
58
59 forAll(toProc, i)
60 {
61 label proci = toProc[i];
62
63 nSend[proci]++;
64 }
65
66
67 // 2. Size sendMap
69
70 forAll(nSend, proci)
71 {
72 sendMap[proci].setSize(nSend[proci]);
73
74 nSend[proci] = 0;
75 }
76
77 // 3. Fill sendMap
78 forAll(toProc, i)
79 {
80 label proci = toProc[i];
81
82 sendMap[proci][nSend[proci]++] = i;
83 }
84
85 // 4. Send over how many I need to receive
86 labelList recvSizes;
87 Pstream::exchangeSizes(sendMap, recvSizes);
88
89
90 // Determine receive map
91 // ~~~~~~~~~~~~~~~~~~~~~
92
93 labelListList constructMap(Pstream::nProcs());
94
95 // Local transfers first
96 constructMap[Pstream::myProcNo()] = identity
97 (
98 sendMap[Pstream::myProcNo()].size()
99 );
100
101 label constructSize = constructMap[Pstream::myProcNo()].size();
102
103 forAll(constructMap, proci)
104 {
105 if (proci != Pstream::myProcNo())
106 {
107 label nRecv = recvSizes[proci];
108
109 constructMap[proci].setSize(nRecv);
110
111 for (label i = 0; i < nRecv; i++)
112 {
113 constructMap[proci][i] = constructSize++;
114 }
115 }
116 }
117
119 (
120 constructSize,
121 std::move(sendMap),
122 std::move(constructMap)
123 );
124}
125
126
127// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
128
129void Foam::backgroundMeshDecomposition::initialRefinement()
130{
131 volScalarField cellWeights
132 (
133 IOobject
134 (
135 "cellWeights",
136 mesh_.time().timeName(),
137 mesh_,
140 ),
141 mesh_,
142 dimensionedScalar("one", dimless, 1.0),
143 zeroGradientFvPatchScalarField::typeName
144 );
145
146 const conformationSurfaces& geometry = geometryToConformTo_;
147
148 decompositionMethod& decomposer =
149 decompositionModel::New(mesh_).decomposer();
150
151 volScalarField::Internal& icellWeights = cellWeights;
152
153 // For each cell in the mesh has it been determined if it is fully
154 // inside, outside, or overlaps the surface
155 List<volumeType> volumeStatus
156 (
157 mesh_.nCells(),
159 );
160
161 // Surface refinement
162 {
163 while (true)
164 {
165 // Determine/update the status of each cell
166 forAll(volumeStatus, celli)
167 {
168 if (volumeStatus[celli] == volumeType::UNKNOWN)
169 {
170 treeBoundBox cellBb
171 (
172 mesh_.cells()[celli].points
173 (
174 mesh_.faces(),
175 mesh_.points()
176 )
177 );
178
179 if (geometry.overlaps(cellBb))
180 {
181 volumeStatus[celli] = volumeType::MIXED;
182 }
183 else if (geometry.inside(cellBb.centre()))
184 {
185 volumeStatus[celli] = volumeType::INSIDE;
186 }
187 else
188 {
189 volumeStatus[celli] = volumeType::OUTSIDE;
190 }
191 }
192 }
193
194 {
195 labelList refCells = selectRefinementCells
196 (
197 volumeStatus,
198 cellWeights
199 );
200
201 // Maintain 2:1 ratio
202 labelList newCellsToRefine
203 (
204 meshCutter_.consistentRefinement
205 (
206 refCells,
207 true // extend set
208 )
209 );
210
211 forAll(newCellsToRefine, nCTRI)
212 {
213 label celli = newCellsToRefine[nCTRI];
214
215 if (volumeStatus[celli] == volumeType::MIXED)
216 {
217 volumeStatus[celli] = volumeType::UNKNOWN;
218 }
219
220 icellWeights[celli] = max
221 (
222 1.0,
223 icellWeights[celli]/8.0
224 );
225 }
226
227 if (returnReduce(newCellsToRefine.size(), sumOp<label>()) == 0)
228 {
229 break;
230 }
231
232 // Mesh changing engine.
233 polyTopoChange meshMod(mesh_);
234
235 // Play refinement commands into mesh changer.
236 meshCutter_.setRefinement(newCellsToRefine, meshMod);
237
238 // Create mesh, return map from old to new mesh.
239 autoPtr<mapPolyMesh> map = meshMod.changeMesh
240 (
241 mesh_,
242 false, // inflate
243 true, // syncParallel
244 true, // orderCells (to reduce cell transfers)
245 false // orderPoints
246 );
247
248 // Update fields
249 mesh_.updateMesh(map());
250
251 // Update numbering of cells/vertices.
252 meshCutter_.updateMesh(map());
253
254 {
255 // Map volumeStatus
256
257 const labelList& cellMap = map().cellMap();
258
259 List<volumeType> newVolumeStatus(cellMap.size());
260
261 forAll(cellMap, newCelli)
262 {
263 label oldCelli = cellMap[newCelli];
264
265 if (oldCelli == -1)
266 {
267 newVolumeStatus[newCelli] = volumeType::UNKNOWN;
268 }
269 else
270 {
271 newVolumeStatus[newCelli] = volumeStatus[oldCelli];
272 }
273 }
274
275 volumeStatus.transfer(newVolumeStatus);
276 }
277
278 Info<< " Background mesh refined from "
279 << returnReduce(map().nOldCells(), sumOp<label>())
280 << " to " << mesh_.globalData().nTotalCells()
281 << " cells." << endl;
282 }
283
284 // Determine/update the status of each cell
285 forAll(volumeStatus, celli)
286 {
287 if (volumeStatus[celli] == volumeType::UNKNOWN)
288 {
289 treeBoundBox cellBb
290 (
291 mesh_.cells()[celli].points
292 (
293 mesh_.faces(),
294 mesh_.points()
295 )
296 );
297
298 if (geometry.overlaps(cellBb))
299 {
300 volumeStatus[celli] = volumeType::MIXED;
301 }
302 else if (geometry.inside(cellBb.centre()))
303 {
304 volumeStatus[celli] = volumeType::INSIDE;
305 }
306 else
307 {
308 volumeStatus[celli] = volumeType::OUTSIDE;
309 }
310 }
311 }
312
313 // Hard code switch for this stage for testing
314 bool removeOutsideCells = false;
315
316 if (removeOutsideCells)
317 {
318 DynamicList<label> cellsToRemove;
319
320 forAll(volumeStatus, celli)
321 {
322 if (volumeStatus[celli] == volumeType::OUTSIDE)
323 {
324 cellsToRemove.append(celli);
325 }
326 }
327
328 removeCells cellRemover(mesh_);
329
330 // Mesh changing engine.
331 polyTopoChange meshMod(mesh_);
332
333 labelList exposedFaces = cellRemover.getExposedFaces
334 (
335 cellsToRemove
336 );
337
338 // Play refinement commands into mesh changer.
339 cellRemover.setRefinement
340 (
341 cellsToRemove,
342 exposedFaces,
343 labelList(exposedFaces.size(), Zero), // patchID dummy
344 meshMod
345 );
346
347 // Create mesh, return map from old to new mesh.
348 autoPtr<mapPolyMesh> map = meshMod.changeMesh
349 (
350 mesh_,
351 false, // inflate
352 true, // syncParallel
353 true, // orderCells (to reduce cell transfers)
354 false // orderPoints
355 );
356
357 // Update fields
358 mesh_.updateMesh(map());
359
360 // Update numbering of cells/vertices.
361 meshCutter_.updateMesh(map());
362 cellRemover.updateMesh(map());
363
364 {
365 // Map volumeStatus
366
367 const labelList& cellMap = map().cellMap();
368
369 List<volumeType> newVolumeStatus(cellMap.size());
370
371 forAll(cellMap, newCelli)
372 {
373 label oldCelli = cellMap[newCelli];
374
375 if (oldCelli == -1)
376 {
377 newVolumeStatus[newCelli] = volumeType::UNKNOWN;
378 }
379 else
380 {
381 newVolumeStatus[newCelli] = volumeStatus[oldCelli];
382 }
383 }
384
385 volumeStatus.transfer(newVolumeStatus);
386 }
387
388 Info<< "Removed "
389 << returnReduce(map().nOldCells(), sumOp<label>())
390 - mesh_.globalData().nTotalCells()
391 << " cells." << endl;
392 }
393
394 if (debug)
395 {
396 // const_cast<Time&>(mesh_.time())++;
397 // Info<< "Time " << mesh_.time().timeName() << endl;
398 meshCutter_.write();
399 mesh_.write();
400 cellWeights.write();
401 }
402
403 labelList newDecomp = decomposer.decompose
404 (
405 mesh_,
406 mesh_.cellCentres(),
407 icellWeights
408 );
409
410 fvMeshDistribute distributor(mesh_);
411
412 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute
413 (
414 newDecomp
415 );
416
417 meshCutter_.distribute(mapDist());
418
419 mapDist().distributeCellData(volumeStatus);
420
421 if (debug)
422 {
423 printMeshData(mesh_);
424
425 // const_cast<Time&>(mesh_.time())++;
426 // Info<< "Time " << mesh_.time().timeName() << endl;
427 meshCutter_.write();
428 mesh_.write();
429 cellWeights.write();
430 }
431 }
432 }
433
434 if (debug)
435 {
436 // const_cast<Time&>(mesh_.time())++;
437 // Info<< "Time " << mesh_.time().timeName() << endl;
438 cellWeights.write();
439 mesh_.write();
440 }
441
442 buildPatchAndTree();
443}
444
445
446void Foam::backgroundMeshDecomposition::printMeshData
447(
448 const polyMesh& mesh
449) const
450{
451 // Collect all data on master
452
453 globalIndex globalCells(mesh.nCells());
454 // labelListList patchNeiProcNo(Pstream::nProcs());
455 // labelListList patchSize(Pstream::nProcs());
456 // const labelList& pPatches = mesh.globalData().processorPatches();
457 // patchNeiProcNo[Pstream::myProcNo()].setSize(pPatches.size());
458 // patchSize[Pstream::myProcNo()].setSize(pPatches.size());
459 // forAll(pPatches, i)
460 // {
461 // const processorPolyPatch& ppp = refCast<const processorPolyPatch>
462 // (
463 // mesh.boundaryMesh()[pPatches[i]]
464 // );
465 // patchNeiProcNo[Pstream::myProcNo()][i] = ppp.neighbProcNo();
466 // patchSize[Pstream::myProcNo()][i] = ppp.size();
467 // }
468 // Pstream::gatherList(patchNeiProcNo);
469 // Pstream::gatherList(patchSize);
470
471
472 // // Print stats
473
474 // globalIndex globalBoundaryFaces(mesh.nBoundaryFaces());
475
476 for (const int proci : Pstream::allProcs())
477 {
478 Info<< "Processor " << proci << " "
479 << "Number of cells = " << globalCells.localSize(proci)
480 << endl;
481
482 // label nProcFaces = 0;
483
484 // const labelList& nei = patchNeiProcNo[proci];
485
486 // forAll(patchNeiProcNo[proci], i)
487 // {
488 // Info<< " Number of faces shared with processor "
489 // << patchNeiProcNo[proci][i] << " = " << patchSize[proci][i]
490 // << endl;
491
492 // nProcFaces += patchSize[proci][i];
493 // }
494
495 // Info<< " Number of processor patches = " << nei.size() << nl
496 // << " Number of processor faces = " << nProcFaces << nl
497 // << " Number of boundary faces = "
498 // << globalBoundaryFaces.localSize(proci) << endl;
499 }
500}
501
502
504(
505 label celli,
506 volumeType volType,
507 scalar& weightEstimate
508) const
509{
510 // Sample the box to find an estimate of the min size, and a volume
511 // estimate when overlapping == true.
512
513// const conformationSurfaces& geometry = geometryToConformTo_;
514
515 treeBoundBox cellBb
516 (
517 mesh_.cells()[celli].points
518 (
519 mesh_.faces(),
520 mesh_.points()
521 )
522 );
523
524 weightEstimate = 1.0;
525
526 if (volType == volumeType::MIXED)
527 {
528// // Assess the cell size at the nearest point on the surface for the
529// // MIXED cells, if the cell is large with respect to the cell size,
530// // then refine it.
531//
532// pointField samplePoints
533// (
534// volRes_*volRes_*volRes_,
535// Zero
536// );
537//
538// // scalar sampleVol = cellBb.volume()/samplePoints.size();
539//
540// vector delta = cellBb.span()/volRes_;
541//
542// label pI = 0;
543//
544// for (label i = 0; i < volRes_; i++)
545// {
546// for (label j = 0; j < volRes_; j++)
547// {
548// for (label k = 0; k < volRes_; k++)
549// {
550// samplePoints[pI++] =
551// cellBb.min()
552// + vector
553// (
554// delta.x()*(i + 0.5),
555// delta.y()*(j + 0.5),
556// delta.z()*(k + 0.5)
557// );
558// }
559// }
560// }
561//
562// List<pointIndexHit> hitInfo;
563// labelList hitSurfaces;
564//
565// geometry.findSurfaceNearest
566// (
567// samplePoints,
568// scalarField(samplePoints.size(), sqr(GREAT)),
569// hitInfo,
570// hitSurfaces
571// );
572//
573// // weightEstimate = 0.0;
574//
575// scalar minCellSize = GREAT;
576//
577// forAll(samplePoints, i)
578// {
579// scalar s = cellShapeControl_.cellSize
580// (
581// hitInfo[i].hitPoint()
582// );
583//
584// // Info<< "cellBb.centre() " << cellBb.centre() << nl
585// // << samplePoints[i] << nl
586// // << hitInfo[i] << nl
587// // << "cellBb.span() " << cellBb.span() << nl
588// // << "cellBb.mag() " << cellBb.mag() << nl
589// // << s << endl;
590//
591// if (s < minCellSize)
592// {
593// minCellSize = max(s, minCellSizeLimit_);
594// }
595//
596// // Estimate the number of points in the cell by the surface size,
597// // this is likely to be too small, so reduce.
598// // weightEstimate += sampleVol/pow3(s);
599// }
600//
601// if (sqr(spanScale_)*sqr(minCellSize) < magSqr(cellBb.span()))
602// {
603// return true;
604// }
605 }
606 else if (volType == volumeType::INSIDE)
607 {
608 // scalar s =
609 // foamyHexMesh_.cellShapeControl_.cellSize(cellBb.centre());
610
611 // Estimate number of points in cell by the size at the cell centre
612 // weightEstimate = cellBb.volume()/pow3(s);
613
614 return false;
615 }
616 // else
617 // {
618 // weightEstimate = 1.0;
619
620 // return false;
621 // }
622
623 return false;
624}
625
626
627Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
628(
629 List<volumeType>& volumeStatus,
630 volScalarField& cellWeights
631) const
632{
633 volScalarField::Internal& icellWeights = cellWeights;
634
635 labelHashSet cellsToRefine;
636
637 // Determine/update the status of each cell
638 forAll(volumeStatus, celli)
639 {
640 if (volumeStatus[celli] == volumeType::MIXED)
641 {
642 if (meshCutter_.cellLevel()[celli] < minLevels_)
643 {
644 cellsToRefine.insert(celli);
645 }
646 }
647
648 if (volumeStatus[celli] != volumeType::OUTSIDE)
649 {
650 if
651 (
652 refineCell
653 (
654 celli,
655 volumeStatus[celli],
656 icellWeights[celli]
657 )
658 )
659 {
660 cellsToRefine.insert(celli);
661 }
662 }
663 }
664
665 return cellsToRefine.toc();
666}
667
668
669void Foam::backgroundMeshDecomposition::buildPatchAndTree()
670{
671 primitivePatch tmpBoundaryFaces
672 (
673 SubList<face>
674 (
675 mesh_.faces(),
676 mesh_.nBoundaryFaces(),
677 mesh_.nInternalFaces()
678 ),
679 mesh_.points()
680 );
681
682 boundaryFacesPtr_.reset
683 (
684 new bPatch
685 (
686 tmpBoundaryFaces.localFaces(),
687 tmpBoundaryFaces.localPoints()
688 )
689 );
690
691 // Overall bb
692 treeBoundBox overallBb(boundaryFacesPtr_().localPoints());
693
694 Random& rnd = rndGen_;
695
696 bFTreePtr_.reset
697 (
698 new indexedOctree<treeDataBPatch>
699 (
701 (
702 false,
703 boundaryFacesPtr_(),
705 ),
706 overallBb.extend(rnd, 1e-4),
707 10, // maxLevel
708 10, // leafSize
709 3.0 // duplicity
710 )
711 );
712
713 // Give the bounds of every processor to every other processor
714 allBackgroundMeshBounds_[Pstream::myProcNo()] = overallBb;
715
716 Pstream::allGatherList(allBackgroundMeshBounds_);
717
718 // find global bounding box
719 globalBackgroundBounds_ = treeBoundBox(boundBox::invertedBox);
720 forAll(allBackgroundMeshBounds_, proci)
721 {
722 globalBackgroundBounds_.add(allBackgroundMeshBounds_[proci]);
723 }
724
725 if (false)
726 {
727 OFstream fStr
728 (
729 mesh_.time().path()
730 /"backgroundMeshDecomposition_proc_"
732 + "_boundaryFaces.obj"
733 );
734
735 const faceList& faces = boundaryFacesPtr_().localFaces();
736 const List<point>& points = boundaryFacesPtr_().localPoints();
737
738 Map<label> foamToObj(points.size());
739
740 label vertI = 0;
741
742 forAll(faces, i)
743 {
744 const face& f = faces[i];
745
746 forAll(f, fPI)
747 {
748 if (foamToObj.insert(f[fPI], vertI))
749 {
750 meshTools::writeOBJ(fStr, points[f[fPI]]);
751 vertI++;
752 }
753 }
754
755 fStr<< 'f';
756
757 forAll(f, fPI)
758 {
759 fStr<< ' ' << foamToObj[f[fPI]] + 1;
760 }
761
762 fStr<< nl;
763 }
764 }
765}
766
767
768// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
769
771(
772 const Time& runTime,
773 Random& rndGen,
774 const conformationSurfaces& geometryToConformTo,
775 const dictionary& coeffsDict,
776 const fileName& decompDictFile
777)
778:
779 runTime_(runTime),
780 geometryToConformTo_(geometryToConformTo),
781 rndGen_(rndGen),
782 mesh_
783 (
784 IOobject
785 (
786 "backgroundMeshDecomposition",
787 runTime_.timeName(),
788 runTime_,
789 IOobject::MUST_READ,
790 IOobject::AUTO_WRITE,
791 false
792 )
793 ),
794 meshCutter_
795 (
796 mesh_,
797 labelList(mesh_.nCells(), Zero),
798 labelList(mesh_.nPoints(), Zero)
799 ),
800 boundaryFacesPtr_(),
801 bFTreePtr_(),
802 allBackgroundMeshBounds_(Pstream::nProcs()),
803 globalBackgroundBounds_(),
804 mergeDist_(1e-6*mesh_.bounds().mag()),
805 spanScale_(coeffsDict.get<scalar>("spanScale")),
806 minCellSizeLimit_
807 (
808 coeffsDict.getOrDefault<scalar>("minCellSizeLimit", 0)
809 ),
810 minLevels_(coeffsDict.get<label>("minLevels")),
811 volRes_(coeffsDict.get<label>("sampleResolution")),
812 maxCellWeightCoeff_(coeffsDict.get<scalar>("maxCellWeightCoeff"))
813{
814 if (!Pstream::parRun())
815 {
817 << "This cannot be used when not running in parallel."
818 << exit(FatalError);
819 }
820
821 const decompositionMethod& decomposer =
822 decompositionModel::New(mesh_, decompDictFile).decomposer();
823
824 if (!decomposer.parallelAware())
825 {
827 << "You have selected decomposition method "
828 << decomposer.typeName
829 << " which is not parallel aware." << endl
830 << exit(FatalError);
831 }
832
833 Info<< nl << "Building initial background mesh decomposition" << endl;
834
835 initialRefinement();
836}
837
838
839// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
840
843(
844 volScalarField& cellWeights
845)
846{
847 if (debug)
848 {
849 // const_cast<Time&>(mesh_.time())++;
850 // Info<< "Time " << mesh_.time().timeName() << endl;
851 cellWeights.write();
852 mesh_.write();
853 }
854
855 volScalarField::Internal& icellWeights = cellWeights;
856
857 while (true)
858 {
859 // Refine large cells if necessary
860
861 label nOccupiedCells = 0;
862
863 forAll(icellWeights, cI)
864 {
865 if (icellWeights[cI] > 1 - SMALL)
866 {
867 nOccupiedCells++;
868 }
869 }
870
871 // Only look at occupied cells, as there is a possibility of runaway
872 // refinement if the number of cells grows too fast. Also, clip the
873 // minimum cellWeightLimit at maxCellWeightCoeff_
874
875 scalar cellWeightLimit = max
876 (
877 maxCellWeightCoeff_
878 *sum(cellWeights).value()
879 /returnReduce(nOccupiedCells, sumOp<label>()),
880 maxCellWeightCoeff_
881 );
882
883 if (debug)
884 {
885 Info<< " cellWeightLimit " << cellWeightLimit << endl;
886
887 Pout<< " sum(cellWeights) " << sum(cellWeights.primitiveField())
888 << " max(cellWeights) " << max(cellWeights.primitiveField())
889 << endl;
890 }
891
892 labelHashSet cellsToRefine;
893
894 forAll(icellWeights, cWI)
895 {
896 if (icellWeights[cWI] > cellWeightLimit)
897 {
898 cellsToRefine.insert(cWI);
899 }
900 }
901
902 if (returnReduce(cellsToRefine.size(), sumOp<label>()) == 0)
903 {
904 break;
905 }
906
907 // Maintain 2:1 ratio
908 labelList newCellsToRefine
909 (
910 meshCutter_.consistentRefinement
911 (
912 cellsToRefine.toc(),
913 true // extend set
914 )
915 );
916
917 if (debug && !cellsToRefine.empty())
918 {
919 Pout<< " cellWeights too large in " << cellsToRefine.size()
920 << " cells" << endl;
921 }
922
923 forAll(newCellsToRefine, nCTRI)
924 {
925 label celli = newCellsToRefine[nCTRI];
926
927 icellWeights[celli] /= 8.0;
928 }
929
930 // Mesh changing engine.
931 polyTopoChange meshMod(mesh_);
932
933 // Play refinement commands into mesh changer.
934 meshCutter_.setRefinement(newCellsToRefine, meshMod);
935
936 // Create mesh, return map from old to new mesh.
937 autoPtr<mapPolyMesh> map = meshMod.changeMesh
938 (
939 mesh_,
940 false, // inflate
941 true, // syncParallel
942 true, // orderCells (to reduce cell motion)
943 false // orderPoints
944 );
945
946 // Update fields
947 mesh_.updateMesh(map());
948
949 // Update numbering of cells/vertices.
950 meshCutter_.updateMesh(map());
951
952 Info<< " Background mesh refined from "
953 << returnReduce(map().nOldCells(), sumOp<label>())
954 << " to " << mesh_.globalData().nTotalCells()
955 << " cells." << endl;
956
957 if (debug)
958 {
959 // const_cast<Time&>(mesh_.time())++;
960 // Info<< "Time " << mesh_.time().timeName() << endl;
961 cellWeights.write();
962 mesh_.write();
963 }
964 }
965
966 if (debug)
967 {
968 printMeshData(mesh_);
969
970 Pout<< " Pre distribute sum(cellWeights) "
971 << sum(icellWeights)
972 << " max(cellWeights) "
973 << max(icellWeights)
974 << endl;
975 }
976
977 decompositionMethod& decomposer =
978 decompositionModel::New(mesh_).decomposer();
979
980 labelList newDecomp = decomposer.decompose
981 (
982 mesh_,
983 mesh_.cellCentres(),
984 icellWeights
985 );
986
987 Info<< " Redistributing background mesh cells" << endl;
988
989 fvMeshDistribute distributor(mesh_);
990
991 autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute(newDecomp);
992
993 meshCutter_.distribute(mapDist());
994
995 if (debug)
996 {
997 printMeshData(mesh_);
998
999 Pout<< " Post distribute sum(cellWeights) "
1000 << sum(icellWeights)
1001 << " max(cellWeights) "
1002 << max(icellWeights)
1003 << endl;
1004
1005 // const_cast<Time&>(mesh_.time())++;
1006 // Info<< "Time " << mesh_.time().timeName() << endl;
1007 mesh_.write();
1008 cellWeights.write();
1009 }
1010
1011 buildPatchAndTree();
1012
1013 return mapDist;
1014}
1015
1016
1018(
1019 const point& pt
1020) const
1021{
1022// return bFTreePtr_().findAnyOverlap(pt, 0.0);
1023
1024 return (bFTreePtr_().getVolumeType(pt) == volumeType::INSIDE);
1025}
1026
1027
1029(
1030 const List<point>& pts
1031) const
1032{
1033 boolList posProc(pts.size(), true);
1034
1035 forAll(pts, pI)
1036 {
1037 posProc[pI] = positionOnThisProcessor(pts[pI]);
1038 }
1039
1040 return posProc;
1041}
1042
1043
1045(
1046 const treeBoundBox& box
1047) const
1048{
1049// return !procBounds().contains(box);
1050 return !bFTreePtr_().findBox(box).empty();
1051}
1052
1053
1055(
1056 const point& centre,
1057 const scalar radiusSqr
1058) const
1059{
1060 //return bFTreePtr_().findAnyOverlap(centre, radiusSqr);
1061
1062 return bFTreePtr_().findNearest(centre, radiusSqr).hit();
1063}
1064
1065
1067(
1068 const point& start,
1069 const point& end
1070) const
1071{
1072 return bFTreePtr_().findLine(start, end);
1073}
1074
1075
1077(
1078 const point& start,
1079 const point& end
1080) const
1081{
1082 return bFTreePtr_().findLineAny(start, end);
1083}
1084
1085
1087(
1088 const List<point>& pts
1089) const
1090{
1091 DynamicList<label> toCandidateProc;
1092 DynamicList<point> testPoints;
1093 labelList ptBlockStart(pts.size(), -1);
1094 labelList ptBlockSize(pts.size(), -1);
1095
1096 label nTotalCandidates = 0;
1097
1098 forAll(pts, pI)
1099 {
1100 const point& pt = pts[pI];
1101
1102 label nCandidates = 0;
1103
1104 forAll(allBackgroundMeshBounds_, proci)
1105 {
1106 // Candidate points may lie just outside a processor box, increase
1107 // test range by using overlaps rather than contains
1108 if (allBackgroundMeshBounds_[proci].overlaps(pt, sqr(SMALL*100)))
1109 {
1110 toCandidateProc.append(proci);
1111 testPoints.append(pt);
1112
1113 nCandidates++;
1114 }
1115 }
1116
1117 ptBlockStart[pI] = nTotalCandidates;
1118 ptBlockSize[pI] = nCandidates;
1119
1120 nTotalCandidates += nCandidates;
1121 }
1122
1123 // Needed for reverseDistribute
1124 label preDistributionToCandidateProcSize = toCandidateProc.size();
1125
1126 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1127
1128 map().distribute(testPoints);
1129
1130 List<scalar> distanceSqrToCandidate(testPoints.size(), sqr(GREAT));
1131
1132 // Test candidate points on candidate processors
1133 forAll(testPoints, tPI)
1134 {
1135 pointIndexHit info = bFTreePtr_().findNearest
1136 (
1137 testPoints[tPI],
1138 sqr(GREAT)
1139 );
1140
1141 if (info.hit())
1142 {
1143 distanceSqrToCandidate[tPI] = magSqr
1144 (
1145 testPoints[tPI] - info.hitPoint()
1146 );
1147 }
1148 }
1149
1150 map().reverseDistribute
1151 (
1152 preDistributionToCandidateProcSize,
1153 distanceSqrToCandidate
1154 );
1155
1156 labelList ptNearestProc(pts.size(), -1);
1157
1158 forAll(pts, pI)
1159 {
1160 // Extract the sub list of results for this point
1161
1162 SubList<scalar> ptNearestProcResults
1163 (
1164 distanceSqrToCandidate,
1165 ptBlockSize[pI],
1166 ptBlockStart[pI]
1167 );
1168
1169 scalar nearestProcDistSqr = GREAT;
1170
1171 forAll(ptNearestProcResults, pPRI)
1172 {
1173 if (ptNearestProcResults[pPRI] < nearestProcDistSqr)
1174 {
1175 nearestProcDistSqr = ptNearestProcResults[pPRI];
1176
1177 ptNearestProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
1178 }
1179 }
1180
1181 if (debug)
1182 {
1183 Pout<< pts[pI] << " nearestProcDistSqr " << nearestProcDistSqr
1184 << " ptNearestProc[pI] " << ptNearestProc[pI] << endl;
1185 }
1186
1187 if (ptNearestProc[pI] < 0)
1188 {
1190 << "The position " << pts[pI]
1191 << " did not find a nearest point on the background mesh."
1192 << exit(FatalError);
1193 }
1194 }
1195
1196 return ptNearestProc;
1197}
1198
1199
1200
1203(
1204 const List<point>& starts,
1205 const List<point>& ends,
1206 bool includeOwnProcessor
1207) const
1208{
1209 DynamicList<label> toCandidateProc;
1210 DynamicList<point> testStarts;
1211 DynamicList<point> testEnds;
1212 labelList segmentBlockStart(starts.size(), -1);
1213 labelList segmentBlockSize(starts.size(), -1);
1214
1215 label nTotalCandidates = 0;
1216
1217 forAll(starts, sI)
1218 {
1219 const point& s = starts[sI];
1220 const point& e = ends[sI];
1221
1222 // Dummy point for treeBoundBox::intersects
1223 point p(Zero);
1224
1225 label nCandidates = 0;
1226
1227 forAll(allBackgroundMeshBounds_, proci)
1228 {
1229 // It is assumed that the sphere in question overlaps the source
1230 // processor, so don't test it, unless includeOwnProcessor is true
1231 if
1232 (
1233 (includeOwnProcessor || proci != Pstream::myProcNo())
1234 && allBackgroundMeshBounds_[proci].intersects(s, e, p)
1235 )
1236 {
1237 toCandidateProc.append(proci);
1238 testStarts.append(s);
1239 testEnds.append(e);
1240
1241 nCandidates++;
1242 }
1243 }
1244
1245 segmentBlockStart[sI] = nTotalCandidates;
1246 segmentBlockSize[sI] = nCandidates;
1247
1248 nTotalCandidates += nCandidates;
1249 }
1250
1251 // Needed for reverseDistribute
1252 label preDistributionToCandidateProcSize = toCandidateProc.size();
1253
1254 autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1255
1256 map().distribute(testStarts);
1257 map().distribute(testEnds);
1258
1259 List<pointIndexHit> segmentIntersectsCandidate(testStarts.size());
1260
1261 // Test candidate segments on candidate processors
1262 forAll(testStarts, sI)
1263 {
1264 const point& s = testStarts[sI];
1265 const point& e = testEnds[sI];
1266
1267 // If the sphere finds some elements of the patch, then it overlaps
1268 segmentIntersectsCandidate[sI] = bFTreePtr_().findLine(s, e);
1269 }
1270
1271 map().reverseDistribute
1272 (
1273 preDistributionToCandidateProcSize,
1274 segmentIntersectsCandidate
1275 );
1276
1277 List<List<pointIndexHit>> segmentHitProcs(starts.size());
1278
1279 // Working storage for assessing processors
1280 DynamicList<pointIndexHit> tmpProcHits;
1281
1282 forAll(starts, sI)
1283 {
1284 tmpProcHits.clear();
1285
1286 // Extract the sub list of results for this point
1287
1288 SubList<pointIndexHit> segmentProcResults
1289 (
1290 segmentIntersectsCandidate,
1291 segmentBlockSize[sI],
1292 segmentBlockStart[sI]
1293 );
1294
1295 forAll(segmentProcResults, sPRI)
1296 {
1297 if (segmentProcResults[sPRI].hit())
1298 {
1299 tmpProcHits.append(segmentProcResults[sPRI]);
1300
1301 tmpProcHits.last().setIndex
1302 (
1303 toCandidateProc[segmentBlockStart[sI] + sPRI]
1304 );
1305 }
1306 }
1307
1308 segmentHitProcs[sI] = tmpProcHits;
1309 }
1310
1311 return segmentHitProcs;
1312}
1313
1314
1316(
1317 const point& centre,
1318 const scalar& radiusSqr
1319) const
1320{
1321 forAll(allBackgroundMeshBounds_, proci)
1322 {
1323 if (bFTreePtr_().findNearest(centre, radiusSqr).hit())
1324 {
1325 return true;
1326 }
1327 }
1328
1329 return false;
1330}
1331
1332
1334(
1335 const point& centre,
1336 const scalar radiusSqr
1337) const
1338{
1339 DynamicList<label> toProc(Pstream::nProcs());
1340
1341 forAll(allBackgroundMeshBounds_, proci)
1342 {
1343 // Test against the bounding box of the processor
1344 if
1345 (
1346 proci != Pstream::myProcNo()
1347 && allBackgroundMeshBounds_[proci].overlaps(centre, radiusSqr)
1348 )
1349 {
1350 // Expensive test
1351// if (bFTreePtr_().findNearest(centre, radiusSqr).hit())
1352 {
1353 toProc.append(proci);
1354 }
1355 }
1356 }
1357
1358 return toProc;
1359}
1360
1361
1362//Foam::labelListList Foam::backgroundMeshDecomposition::overlapsProcessors
1363//(
1364// const List<point>& centres,
1365// const List<scalar>& radiusSqrs,
1366// const Delaunay& T,
1367// bool includeOwnProcessor
1368//) const
1369//{
1370// DynamicList<label> toCandidateProc;
1371// DynamicList<point> testCentres;
1372// DynamicList<scalar> testRadiusSqrs;
1373// labelList sphereBlockStart(centres.size(), -1);
1374// labelList sphereBlockSize(centres.size(), -1);
1375//
1376// label nTotalCandidates = 0;
1377//
1378// forAll(centres, sI)
1379// {
1380// const point& c = centres[sI];
1381// scalar rSqr = radiusSqrs[sI];
1382//
1383// label nCandidates = 0;
1384//
1385// forAll(allBackgroundMeshBounds_, proci)
1386// {
1387// // It is assumed that the sphere in question overlaps the source
1388// // processor, so don't test it, unless includeOwnProcessor is true
1389// if
1390// (
1391// (includeOwnProcessor || proci != Pstream::myProcNo())
1392// && allBackgroundMeshBounds_[proci].overlaps(c, rSqr)
1393// )
1394// {
1395// if (bFTreePtr_().findNearest(c, rSqr).hit())
1396// {
1397// toCandidateProc.append(proci);
1398// testCentres.append(c);
1399// testRadiusSqrs.append(rSqr);
1400//
1401// nCandidates++;
1402// }
1403// }
1404// }
1405//
1406// sphereBlockStart[sI] = nTotalCandidates;
1407// sphereBlockSize[sI] = nCandidates;
1408//
1409// nTotalCandidates += nCandidates;
1410// }
1411//
1412// // Needed for reverseDistribute
1419//
1420// // TODO: This is faster, but results in more vertices being referred
1421// boolList sphereOverlapsCandidate(testCentres.size(), true);
1458//
1464//
1465// labelListList sphereProcs(centres.size());
1466//
1467// // Working storage for assessing processors
1468// DynamicList<label> tmpProcs;
1469//
1470// forAll(centres, sI)
1471// {
1472// tmpProcs.clear();
1473//
1474// // Extract the sub list of results for this point
1475//
1476// SubList<bool> sphereProcResults
1477// (
1478// sphereOverlapsCandidate,
1479// sphereBlockSize[sI],
1480// sphereBlockStart[sI]
1481// );
1482//
1483// forAll(sphereProcResults, sPRI)
1484// {
1485// if (sphereProcResults[sPRI])
1486// {
1487// tmpProcs.append(toCandidateProc[sphereBlockStart[sI] + sPRI]);
1488// }
1489// }
1490//
1491// sphereProcs[sI] = tmpProcs;
1492// }
1493//
1494// return sphereProcs;
1495//}
1496
1497
1498//Foam::labelListList Foam::backgroundMeshDecomposition::overlapProcessors
1499//(
1500// const point& cc,
1501// const scalar rSqr
1502//) const
1503//{
1504// DynamicList<label> toCandidateProc;
1505// label sphereBlockStart(-1);
1506// label sphereBlockSize(-1);
1507//
1508// label nCandidates = 0;
1509//
1510// forAll(allBackgroundMeshBounds_, proci)
1511// {
1512// // It is assumed that the sphere in question overlaps the source
1513// // processor, so don't test it, unless includeOwnProcessor is true
1514// if
1515// (
1516// (includeOwnProcessor || proci != Pstream::myProcNo())
1517// && allBackgroundMeshBounds_[proci].overlaps(cc, rSqr)
1518// )
1519// {
1520// toCandidateProc.append(proci);
1521//
1522// nCandidates++;
1523// }
1524// }
1525//
1526// sphereBlockSize = nCandidates;
1527// nTotalCandidates += nCandidates;
1528//
1529// // Needed for reverseDistribute
1530// label preDistributionToCandidateProcSize = toCandidateProc.size();
1531//
1532// autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1533//
1534// map().distribute(testCentres);
1535// map().distribute(testRadiusSqrs);
1536//
1537// // TODO: This is faster, but results in more vertices being referred
1539// boolList sphereOverlapsCandidate(testCentres.size(), false);
1540//
1541// // Test candidate spheres on candidate processors
1542// forAll(testCentres, sI)
1543// {
1544// const point& c = testCentres[sI];
1545// const scalar rSqr = testRadiusSqrs[sI];
1546//
1547// const bool flagOverlap = bFTreePtr_().findNearest(c, rSqr).hit();
1548//
1549// if (flagOverlap)
1550// {
1551// //if (vertexOctree.findAnyOverlap(c, rSqr))
1556//
1571// sphereOverlapsCandidate[sI] = true;
1573// }
1574// }
1575//
1576// map().reverseDistribute
1577// (
1578// preDistributionToCandidateProcSize,
1579// sphereOverlapsCandidate
1580// );
1581//
1582// labelListList sphereProcs(centres.size());
1583//
1584// // Working storage for assessing processors
1585// DynamicList<label> tmpProcs;
1586//
1587// forAll(centres, sI)
1588// {
1589// tmpProcs.clear();
1590//
1591// // Extract the sub list of results for this point
1592//
1593// SubList<bool> sphereProcResults
1594// (
1595// sphereOverlapsCandidate,
1596// sphereBlockSize[sI],
1597// sphereBlockStart[sI]
1598// );
1599//
1600// forAll(sphereProcResults, sPRI)
1601// {
1602// if (sphereProcResults[sPRI])
1603// {
1604// tmpProcs.append(toCandidateProc[sphereBlockStart[sI] + sPRI]);
1605// }
1606// }
1607//
1608// sphereProcs[sI] = tmpProcs;
1609// }
1610//
1611// return sphereProcs;
1612//}
1613
1614
1615// ************************************************************************* //
DimensionedField< scalar, volMesh > Internal
The internal field type from which this GeometricField is derived.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
UPstream::rangeType allProcs() const noexcept
Range of ranks indices associated with PstreamBuffers.
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 exchangeSizes(const labelUList &sendProcs, const labelUList &recvProcs, const Container &sendData, labelList &sizes, const label tag=UPstream::msgType(), const label comm=UPstream::worldComm)
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
Store a background polyMesh to use for the decomposition of space and queries for parallel conformalV...
labelList processorNearestPosition(const List< point > &pts) const
What is the nearest processor to the given position?
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end, (exposing.
bool overlapsThisProcessor(const treeBoundBox &box) const
Does the given box overlap the faces of the boundary of this.
bool positionOnThisProcessor(const point &pt) const
Is the given position inside the domain of this decomposition.
autoPtr< mapDistributePolyMesh > distribute(volScalarField &cellWeights)
Redistribute the background mesh based on a supplied weight field,.
static autoPtr< mapDistribute > buildMap(const List< label > &toProc)
Build a mapDistribute for the supplied destination processor data.
List< List< pointIndexHit > > intersectsProcessors(const List< point > &starts, const List< point > &ends, bool includeOwnProcessor=false) const
Which processors are intersected by the line segment, returns all.
pointIndexHit findLine(const point &start, const point &end) const
Find nearest intersection of line between start and end, (exposing.
bool overlapsOtherProcessors(const point &centre, const scalar &radiusSqr) const
labelList overlapProcessors(const point &centre, const scalar radiusSqr) const
static const boundBox invertedBox
A large inverted boundBox: min/max == +/- ROOTVGREAT.
Definition: boundBox.H:86
virtual bool write(const bool valid=true) const
Write mesh using IO settings from time.
Definition: fvMesh.C:1079
virtual void updateMesh(const mapPolyMesh &mpm)
Update mesh corresponding to the given map.
Definition: fvMesh.C:971
label nTotalCells() const noexcept
Return total number of cells in decomposed mesh.
labelListList setRefinement(const labelList &cells, polyTopoChange &)
Insert refinement. All selected cells will be split into 8.
Definition: hexRef8.C:3190
labelList consistentRefinement(const labelUList &cellLevel, const labelList &cellsToRefine, const bool maxSet) const
Given valid mesh and current cell level and proposed.
Definition: hexRef8.C:2250
void distribute(const mapDistributePolyMesh &)
Update local numbering for mesh redistribution.
Definition: hexRef8.C:4516
void updateMesh(const mapPolyMesh &)
Update local numbering for changed mesh.
Definition: hexRef8.C:4231
static scalar & perturbTol()
Get the perturbation tolerance.
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1310
const vectorField & cellCentres() const
label nCells() const noexcept
Number of mesh cells.
int myProcNo() const noexcept
Return processor number.
Container with cells to refine. Refinement given as single direction.
Definition: refineCell.H:57
@ OUTSIDE
A location outside the volume.
Definition: volumeType.H:69
@ MIXED
A location that is partly inside and outside.
Definition: volumeType.H:70
@ UNKNOWN
Unknown state.
Definition: volumeType.H:67
@ INSIDE
A location inside the volume.
Definition: volumeType.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
volScalarField & p
dynamicFvMesh & mesh
engineTime & runTime
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
label nPoints
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
word timeName
Definition: getTimeIndex.H:3
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
Namespace for OpenFOAM.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
const dimensionSet dimless
Dimensionless.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
List< label > labelList
A List of labels.
Definition: List.H:66
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
vector point
Point is a vector.
Definition: point.H:43
treeDataPrimitivePatch< bPatch > treeDataBPatch
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
List< bool > boolList
A List of bools.
Definition: List.H:64
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
PrimitivePatch<::Foam::List< face >, const pointField > bPatch
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
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)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Random rndGen
Definition: createFields.H:23