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