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 | Copyright (C) 2011-2016 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"
30 #include "Time.H"
31 #include "Random.H"
32 #include "pointConversion.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 
39 defineTypeNameAndDebug(backgroundMeshDecomposition, 0);
40 
41 }
42 
43 
44 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
45 
47 (
48  const List<label>& toProc
49 )
50 {
51  // Determine send map
52  // ~~~~~~~~~~~~~~~~~~
53 
54  // 1. Count
55  labelList nSend(Pstream::nProcs(), 0);
56 
57  forAll(toProc, i)
58  {
59  label proci = toProc[i];
60 
61  nSend[proci]++;
62  }
63 
64 
65  // 2. Size sendMap
66  labelListList sendMap(Pstream::nProcs());
67 
68  forAll(nSend, proci)
69  {
70  sendMap[proci].setSize(nSend[proci]);
71 
72  nSend[proci] = 0;
73  }
74 
75  // 3. Fill sendMap
76  forAll(toProc, i)
77  {
78  label proci = toProc[i];
79 
80  sendMap[proci][nSend[proci]++] = i;
81  }
82 
83  // 4. Send over how many I need to receive
84  labelList recvSizes;
85  Pstream::exchangeSizes(sendMap, recvSizes);
86 
87 
88  // Determine receive map
89  // ~~~~~~~~~~~~~~~~~~~~~
90 
91  labelListList constructMap(Pstream::nProcs());
92 
93  // Local transfers first
94  constructMap[Pstream::myProcNo()] = identity
95  (
96  sendMap[Pstream::myProcNo()].size()
97  );
98 
99  label constructSize = constructMap[Pstream::myProcNo()].size();
100 
101  forAll(constructMap, proci)
102  {
103  if (proci != Pstream::myProcNo())
104  {
105  label nRecv = recvSizes[proci];
106 
107  constructMap[proci].setSize(nRecv);
108 
109  for (label i = 0; i < nRecv; i++)
110  {
111  constructMap[proci][i] = constructSize++;
112  }
113  }
114  }
115 
116  return autoPtr<mapDistribute>
117  (
118  new mapDistribute
119  (
120  constructSize,
121  sendMap.xfer(),
122  constructMap.xfer()
123  )
124  );
125 }
126 
127 
128 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
129 
130 void Foam::backgroundMeshDecomposition::initialRefinement()
131 {
132  volScalarField cellWeights
133  (
134  IOobject
135  (
136  "cellWeights",
137  mesh_.time().timeName(),
138  mesh_,
141  ),
142  mesh_,
143  dimensionedScalar("one", dimless, 1.0),
144  zeroGradientFvPatchScalarField::typeName
145  );
146 
147  const conformationSurfaces& geometry = geometryToConformTo_;
148 
149  decompositionMethod& decomposer = decomposerPtr_();
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.midpoint()))
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.midpoint()))
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(), 0), // 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_, mergeDist_);
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 
446 void 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.nFaces()-mesh.nInternalFaces());
475 
476  for (label proci = 0; proci < Pstream::nProcs(); proci++)
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 
503 bool Foam::backgroundMeshDecomposition::refineCell
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.midpoint() " << cellBb.midpoint() << 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.midpoint());
610 
611  // Estimate the number of points in the cell by the size at the cell
612  // midpoint
613  // weightEstimate = cellBb.volume()/pow3(s);
614 
615  return false;
616  }
617  // else
618  // {
619  // weightEstimate = 1.0;
620 
621  // return false;
622  // }
623 
624  return false;
625 }
626 
627 
628 Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
629 (
630  List<volumeType>& volumeStatus,
631  volScalarField& cellWeights
632 ) const
633 {
634  volScalarField::Internal& icellWeights = cellWeights;
635 
636  labelHashSet cellsToRefine;
637 
638  // Determine/update the status of each cell
639  forAll(volumeStatus, celli)
640  {
641  if (volumeStatus[celli] == volumeType::MIXED)
642  {
643  if (meshCutter_.cellLevel()[celli] < minLevels_)
644  {
645  cellsToRefine.insert(celli);
646  }
647  }
648 
649  if (volumeStatus[celli] != volumeType::OUTSIDE)
650  {
651  if
652  (
653  refineCell
654  (
655  celli,
656  volumeStatus[celli],
657  icellWeights[celli]
658  )
659  )
660  {
661  cellsToRefine.insert(celli);
662  }
663  }
664  }
665 
666  return cellsToRefine.toc();
667 }
668 
669 
670 void Foam::backgroundMeshDecomposition::buildPatchAndTree()
671 {
672  primitivePatch tmpBoundaryFaces
673  (
674  SubList<face>
675  (
676  mesh_.faces(),
677  mesh_.nFaces() - mesh_.nInternalFaces(),
678  mesh_.nInternalFaces()
679  ),
680  mesh_.points()
681  );
682 
683  boundaryFacesPtr_.reset
684  (
685  new bPatch
686  (
687  tmpBoundaryFaces.localFaces(),
688  tmpBoundaryFaces.localPoints()
689  )
690  );
691 
692  // Overall bb
693  treeBoundBox overallBb(boundaryFacesPtr_().localPoints());
694 
695  Random& rnd = rndGen_;
696 
697  bFTreePtr_.reset
698  (
699  new indexedOctree<treeDataBPatch>
700  (
702  (
703  false,
704  boundaryFacesPtr_(),
706  ),
707  overallBb.extend(rnd, 1e-4),
708  10, // maxLevel
709  10, // leafSize
710  3.0 // duplicity
711  )
712  );
713 
714  // Give the bounds of every processor to every other processor
715  allBackgroundMeshBounds_[Pstream::myProcNo()] = overallBb;
716 
717  Pstream::gatherList(allBackgroundMeshBounds_);
718  Pstream::scatterList(allBackgroundMeshBounds_);
719 
720  point bbMin(GREAT, GREAT, GREAT);
721  point bbMax(-GREAT, -GREAT, -GREAT);
722 
723  forAll(allBackgroundMeshBounds_, proci)
724  {
725  bbMin = min(bbMin, allBackgroundMeshBounds_[proci].min());
726  bbMax = max(bbMax, allBackgroundMeshBounds_[proci].max());
727  }
728 
729  globalBackgroundBounds_ = treeBoundBox(bbMin, bbMax);
730 
731  if (false)
732  {
733  OFstream fStr
734  (
735  mesh_.time().path()
736  /"backgroundMeshDecomposition_proc_"
738  + "_boundaryFaces.obj"
739  );
740 
741  const faceList& faces = boundaryFacesPtr_().localFaces();
742  const List<point>& points = boundaryFacesPtr_().localPoints();
743 
744  Map<label> foamToObj(points.size());
745 
746  label vertI = 0;
747 
748  forAll(faces, i)
749  {
750  const face& f = faces[i];
751 
752  forAll(f, fPI)
753  {
754  if (foamToObj.insert(f[fPI], vertI))
755  {
756  meshTools::writeOBJ(fStr, points[f[fPI]]);
757  vertI++;
758  }
759  }
760 
761  fStr<< 'f';
762 
763  forAll(f, fPI)
764  {
765  fStr<< ' ' << foamToObj[f[fPI]] + 1;
766  }
767 
768  fStr<< nl;
769  }
770  }
771 }
772 
773 
774 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
775 
776 Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
777 (
778  const Time& runTime,
779  Random& rndGen,
780  const conformationSurfaces& geometryToConformTo,
781  const dictionary& coeffsDict
782 )
783 :
784  runTime_(runTime),
785  geometryToConformTo_(geometryToConformTo),
786  rndGen_(rndGen),
787  mesh_
788  (
789  IOobject
790  (
791  "backgroundMeshDecomposition",
792  runTime_.timeName(),
793  runTime_,
794  IOobject::MUST_READ,
795  IOobject::AUTO_WRITE,
796  false
797  )
798  ),
799  meshCutter_
800  (
801  mesh_,
802  labelList(mesh_.nCells(), 0),
803  labelList(mesh_.nPoints(), 0)
804  ),
805  boundaryFacesPtr_(),
806  bFTreePtr_(),
807  allBackgroundMeshBounds_(Pstream::nProcs()),
808  globalBackgroundBounds_(),
809  decomposeDict_
810  (
811  IOobject
812  (
813  "decomposeParDict",
814  runTime_.system(),
815  runTime_,
816  IOobject::MUST_READ_IF_MODIFIED,
817  IOobject::NO_WRITE
818  )
819  ),
820  decomposerPtr_(decompositionMethod::New(decomposeDict_)),
821  mergeDist_(1e-6*mesh_.bounds().mag()),
822  spanScale_(readScalar(coeffsDict.lookup("spanScale"))),
823  minCellSizeLimit_
824  (
825  coeffsDict.lookupOrDefault<scalar>("minCellSizeLimit", 0.0)
826  ),
827  minLevels_(readLabel(coeffsDict.lookup("minLevels"))),
828  volRes_(readLabel(coeffsDict.lookup("sampleResolution"))),
829  maxCellWeightCoeff_(readScalar(coeffsDict.lookup("maxCellWeightCoeff")))
830 {
831  if (!Pstream::parRun())
832  {
834  << "This cannot be used when not running in parallel."
835  << exit(FatalError);
836  }
837 
838  if (!decomposerPtr_().parallelAware())
839  {
841  << "You have selected decomposition method "
842  << decomposerPtr_().typeName
843  << " which is not parallel aware." << endl
844  << exit(FatalError);
845  }
846 
847  Info<< nl << "Building initial background mesh decomposition" << endl;
848 
849  initialRefinement();
850 }
851 
852 
853 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
854 
856 {}
857 
858 
859 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
860 
863 (
864  volScalarField& cellWeights
865 )
866 {
867  if (debug)
868  {
869  // const_cast<Time&>(mesh_.time())++;
870  // Info<< "Time " << mesh_.time().timeName() << endl;
871  cellWeights.write();
872  mesh_.write();
873  }
874 
875  volScalarField::Internal& icellWeights = cellWeights;
876 
877  while (true)
878  {
879  // Refine large cells if necessary
880 
881  label nOccupiedCells = 0;
882 
883  forAll(icellWeights, cI)
884  {
885  if (icellWeights[cI] > 1 - SMALL)
886  {
887  nOccupiedCells++;
888  }
889  }
890 
891  // Only look at occupied cells, as there is a possibility of runaway
892  // refinement if the number of cells grows too fast. Also, clip the
893  // minimum cellWeightLimit at maxCellWeightCoeff_
894 
895  scalar cellWeightLimit = max
896  (
897  maxCellWeightCoeff_
898  *sum(cellWeights).value()
899  /returnReduce(nOccupiedCells, sumOp<label>()),
900  maxCellWeightCoeff_
901  );
902 
903  if (debug)
904  {
905  Info<< " cellWeightLimit " << cellWeightLimit << endl;
906 
907  Pout<< " sum(cellWeights) " << sum(cellWeights.primitiveField())
908  << " max(cellWeights) " << max(cellWeights.primitiveField())
909  << endl;
910  }
911 
912  labelHashSet cellsToRefine;
913 
914  forAll(icellWeights, cWI)
915  {
916  if (icellWeights[cWI] > cellWeightLimit)
917  {
918  cellsToRefine.insert(cWI);
919  }
920  }
921 
922  if (returnReduce(cellsToRefine.size(), sumOp<label>()) == 0)
923  {
924  break;
925  }
926 
927  // Maintain 2:1 ratio
928  labelList newCellsToRefine
929  (
930  meshCutter_.consistentRefinement
931  (
932  cellsToRefine.toc(),
933  true // extend set
934  )
935  );
936 
937  if (debug && !cellsToRefine.empty())
938  {
939  Pout<< " cellWeights too large in " << cellsToRefine.size()
940  << " cells" << endl;
941  }
942 
943  forAll(newCellsToRefine, nCTRI)
944  {
945  label celli = newCellsToRefine[nCTRI];
946 
947  icellWeights[celli] /= 8.0;
948  }
949 
950  // Mesh changing engine.
951  polyTopoChange meshMod(mesh_);
952 
953  // Play refinement commands into mesh changer.
954  meshCutter_.setRefinement(newCellsToRefine, meshMod);
955 
956  // Create mesh, return map from old to new mesh.
957  autoPtr<mapPolyMesh> map = meshMod.changeMesh
958  (
959  mesh_,
960  false, // inflate
961  true, // syncParallel
962  true, // orderCells (to reduce cell motion)
963  false // orderPoints
964  );
965 
966  // Update fields
967  mesh_.updateMesh(map);
968 
969  // Update numbering of cells/vertices.
970  meshCutter_.updateMesh(map);
971 
972  Info<< " Background mesh refined from "
973  << returnReduce(map().nOldCells(), sumOp<label>())
974  << " to " << mesh_.globalData().nTotalCells()
975  << " cells." << endl;
976 
977  if (debug)
978  {
979  // const_cast<Time&>(mesh_.time())++;
980  // Info<< "Time " << mesh_.time().timeName() << endl;
981  cellWeights.write();
982  mesh_.write();
983  }
984  }
985 
986  if (debug)
987  {
988  printMeshData(mesh_);
989 
990  Pout<< " Pre distribute sum(cellWeights) "
991  << sum(icellWeights)
992  << " max(cellWeights) "
993  << max(icellWeights)
994  << endl;
995  }
996 
997  labelList newDecomp = decomposerPtr_().decompose
998  (
999  mesh_,
1000  mesh_.cellCentres(),
1001  icellWeights
1002  );
1003 
1004  Info<< " Redistributing background mesh cells" << endl;
1005 
1006  fvMeshDistribute distributor(mesh_, mergeDist_);
1007 
1008  autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute(newDecomp);
1009 
1010  meshCutter_.distribute(mapDist);
1011 
1012  if (debug)
1013  {
1014  printMeshData(mesh_);
1015 
1016  Pout<< " Post distribute sum(cellWeights) "
1017  << sum(icellWeights)
1018  << " max(cellWeights) "
1019  << max(icellWeights)
1020  << endl;
1021 
1022  // const_cast<Time&>(mesh_.time())++;
1023  // Info<< "Time " << mesh_.time().timeName() << endl;
1024  mesh_.write();
1025  cellWeights.write();
1026  }
1027 
1028  buildPatchAndTree();
1029 
1030  return mapDist;
1031 }
1032 
1033 
1035 (
1036  const point& pt
1037 ) const
1038 {
1039 // return bFTreePtr_().findAnyOverlap(pt, 0.0);
1040 
1041  return (bFTreePtr_().getVolumeType(pt) == volumeType::INSIDE);
1042 }
1043 
1044 
1046 (
1047  const List<point>& pts
1048 ) const
1049 {
1050  boolList posProc(pts.size(), true);
1051 
1052  forAll(pts, pI)
1053  {
1054  posProc[pI] = positionOnThisProcessor(pts[pI]);
1055  }
1056 
1057  return posProc;
1058 }
1059 
1060 
1062 (
1063  const treeBoundBox& box
1064 ) const
1065 {
1066 // return !procBounds().contains(box);
1067  return !bFTreePtr_().findBox(box).empty();
1068 }
1069 
1070 
1072 (
1073  const point& centre,
1074  const scalar radiusSqr
1075 ) const
1076 {
1077  //return bFTreePtr_().findAnyOverlap(centre, radiusSqr);
1078 
1079  return bFTreePtr_().findNearest(centre, radiusSqr).hit();
1080 }
1081 
1082 
1084 (
1085  const point& start,
1086  const point& end
1087 ) const
1088 {
1089  return bFTreePtr_().findLine(start, end);
1090 }
1091 
1092 
1094 (
1095  const point& start,
1096  const point& end
1097 ) const
1098 {
1099  return bFTreePtr_().findLineAny(start, end);
1100 }
1101 
1102 
1104 (
1105  const List<point>& pts
1106 ) const
1107 {
1108  DynamicList<label> toCandidateProc;
1109  DynamicList<point> testPoints;
1110  labelList ptBlockStart(pts.size(), -1);
1111  labelList ptBlockSize(pts.size(), -1);
1112 
1113  label nTotalCandidates = 0;
1114 
1115  forAll(pts, pI)
1116  {
1117  const point& pt = pts[pI];
1118 
1119  label nCandidates = 0;
1120 
1121  forAll(allBackgroundMeshBounds_, proci)
1122  {
1123  // Candidate points may lie just outside a processor box, increase
1124  // test range by using overlaps rather than contains
1125  if (allBackgroundMeshBounds_[proci].overlaps(pt, sqr(SMALL*100)))
1126  {
1127  toCandidateProc.append(proci);
1128  testPoints.append(pt);
1129 
1130  nCandidates++;
1131  }
1132  }
1133 
1134  ptBlockStart[pI] = nTotalCandidates;
1135  ptBlockSize[pI] = nCandidates;
1136 
1137  nTotalCandidates += nCandidates;
1138  }
1139 
1140  // Needed for reverseDistribute
1141  label preDistributionToCandidateProcSize = toCandidateProc.size();
1142 
1143  autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1144 
1145  map().distribute(testPoints);
1146 
1147  List<scalar> distanceSqrToCandidate(testPoints.size(), sqr(GREAT));
1148 
1149  // Test candidate points on candidate processors
1150  forAll(testPoints, tPI)
1151  {
1152  pointIndexHit info = bFTreePtr_().findNearest
1153  (
1154  testPoints[tPI],
1155  sqr(GREAT)
1156  );
1157 
1158  if (info.hit())
1159  {
1160  distanceSqrToCandidate[tPI] = magSqr
1161  (
1162  testPoints[tPI] - info.hitPoint()
1163  );
1164  }
1165  }
1166 
1167  map().reverseDistribute
1168  (
1169  preDistributionToCandidateProcSize,
1170  distanceSqrToCandidate
1171  );
1172 
1173  labelList ptNearestProc(pts.size(), -1);
1174 
1175  forAll(pts, pI)
1176  {
1177  // Extract the sub list of results for this point
1178 
1179  SubList<scalar> ptNearestProcResults
1180  (
1181  distanceSqrToCandidate,
1182  ptBlockSize[pI],
1183  ptBlockStart[pI]
1184  );
1185 
1186  scalar nearestProcDistSqr = GREAT;
1187 
1188  forAll(ptNearestProcResults, pPRI)
1189  {
1190  if (ptNearestProcResults[pPRI] < nearestProcDistSqr)
1191  {
1192  nearestProcDistSqr = ptNearestProcResults[pPRI];
1193 
1194  ptNearestProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
1195  }
1196  }
1197 
1198  if (debug)
1199  {
1200  Pout<< pts[pI] << " nearestProcDistSqr " << nearestProcDistSqr
1201  << " ptNearestProc[pI] " << ptNearestProc[pI] << endl;
1202  }
1203 
1204  if (ptNearestProc[pI] < 0)
1205  {
1207  << "The position " << pts[pI]
1208  << " did not find a nearest point on the background mesh."
1209  << exit(FatalError);
1210  }
1211  }
1212 
1213  return ptNearestProc;
1214 }
1215 
1216 
1217 
1220 (
1221  const List<point>& starts,
1222  const List<point>& ends,
1223  bool includeOwnProcessor
1224 ) const
1225 {
1226  DynamicList<label> toCandidateProc;
1227  DynamicList<point> testStarts;
1228  DynamicList<point> testEnds;
1229  labelList segmentBlockStart(starts.size(), -1);
1230  labelList segmentBlockSize(starts.size(), -1);
1231 
1232  label nTotalCandidates = 0;
1233 
1234  forAll(starts, sI)
1235  {
1236  const point& s = starts[sI];
1237  const point& e = ends[sI];
1238 
1239  // Dummy point for treeBoundBox::intersects
1240  point p(Zero);
1241 
1242  label nCandidates = 0;
1243 
1244  forAll(allBackgroundMeshBounds_, proci)
1245  {
1246  // It is assumed that the sphere in question overlaps the source
1247  // processor, so don't test it, unless includeOwnProcessor is true
1248  if
1249  (
1250  (includeOwnProcessor || proci != Pstream::myProcNo())
1251  && allBackgroundMeshBounds_[proci].intersects(s, e, p)
1252  )
1253  {
1254  toCandidateProc.append(proci);
1255  testStarts.append(s);
1256  testEnds.append(e);
1257 
1258  nCandidates++;
1259  }
1260  }
1261 
1262  segmentBlockStart[sI] = nTotalCandidates;
1263  segmentBlockSize[sI] = nCandidates;
1264 
1265  nTotalCandidates += nCandidates;
1266  }
1267 
1268  // Needed for reverseDistribute
1269  label preDistributionToCandidateProcSize = toCandidateProc.size();
1270 
1271  autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1272 
1273  map().distribute(testStarts);
1274  map().distribute(testEnds);
1275 
1276  List<pointIndexHit> segmentIntersectsCandidate(testStarts.size());
1277 
1278  // Test candidate segments on candidate processors
1279  forAll(testStarts, sI)
1280  {
1281  const point& s = testStarts[sI];
1282  const point& e = testEnds[sI];
1283 
1284  // If the sphere finds some elements of the patch, then it overlaps
1285  segmentIntersectsCandidate[sI] = bFTreePtr_().findLine(s, e);
1286  }
1287 
1288  map().reverseDistribute
1289  (
1290  preDistributionToCandidateProcSize,
1291  segmentIntersectsCandidate
1292  );
1293 
1294  List<List<pointIndexHit>> segmentHitProcs(starts.size());
1295 
1296  // Working storage for assessing processors
1297  DynamicList<pointIndexHit> tmpProcHits;
1298 
1299  forAll(starts, sI)
1300  {
1301  tmpProcHits.clear();
1302 
1303  // Extract the sub list of results for this point
1304 
1305  SubList<pointIndexHit> segmentProcResults
1306  (
1307  segmentIntersectsCandidate,
1308  segmentBlockSize[sI],
1309  segmentBlockStart[sI]
1310  );
1311 
1312  forAll(segmentProcResults, sPRI)
1313  {
1314  if (segmentProcResults[sPRI].hit())
1315  {
1316  tmpProcHits.append(segmentProcResults[sPRI]);
1317 
1318  tmpProcHits.last().setIndex
1319  (
1320  toCandidateProc[segmentBlockStart[sI] + sPRI]
1321  );
1322  }
1323  }
1324 
1325  segmentHitProcs[sI] = tmpProcHits;
1326  }
1327 
1328  return segmentHitProcs;
1329 }
1330 
1331 
1333 (
1334  const point& centre,
1335  const scalar& radiusSqr
1336 ) const
1337 {
1338  forAll(allBackgroundMeshBounds_, proci)
1339  {
1340  if (bFTreePtr_().findNearest(centre, radiusSqr).hit())
1341  {
1342  return true;
1343  }
1344  }
1345 
1346  return false;
1347 }
1348 
1349 
1351 (
1352  const point& centre,
1353  const scalar radiusSqr
1354 ) const
1355 {
1356  DynamicList<label> toProc(Pstream::nProcs());
1357 
1358  forAll(allBackgroundMeshBounds_, proci)
1359  {
1360  // Test against the bounding box of the processor
1361  if
1362  (
1363  proci != Pstream::myProcNo()
1364  && allBackgroundMeshBounds_[proci].overlaps(centre, radiusSqr)
1365  )
1366  {
1367  // Expensive test
1368 // if (bFTreePtr_().findNearest(centre, radiusSqr).hit())
1369  {
1370  toProc.append(proci);
1371  }
1372  }
1373  }
1374 
1375  return toProc;
1376 }
1377 
1378 
1379 //Foam::labelListList Foam::backgroundMeshDecomposition::overlapsProcessors
1380 //(
1381 // const List<point>& centres,
1382 // const List<scalar>& radiusSqrs,
1383 // const Delaunay& T,
1384 // bool includeOwnProcessor
1385 //) const
1386 //{
1387 // DynamicList<label> toCandidateProc;
1388 // DynamicList<point> testCentres;
1389 // DynamicList<scalar> testRadiusSqrs;
1390 // labelList sphereBlockStart(centres.size(), -1);
1391 // labelList sphereBlockSize(centres.size(), -1);
1392 //
1393 // label nTotalCandidates = 0;
1394 //
1395 // forAll(centres, sI)
1396 // {
1397 // const point& c = centres[sI];
1398 // scalar rSqr = radiusSqrs[sI];
1399 //
1400 // label nCandidates = 0;
1401 //
1402 // forAll(allBackgroundMeshBounds_, proci)
1403 // {
1404 // // It is assumed that the sphere in question overlaps the source
1405 // // processor, so don't test it, unless includeOwnProcessor is true
1406 // if
1407 // (
1408 // (includeOwnProcessor || proci != Pstream::myProcNo())
1409 // && allBackgroundMeshBounds_[proci].overlaps(c, rSqr)
1410 // )
1411 // {
1412 // if (bFTreePtr_().findNearest(c, rSqr).hit())
1413 // {
1414 // toCandidateProc.append(proci);
1415 // testCentres.append(c);
1416 // testRadiusSqrs.append(rSqr);
1417 //
1418 // nCandidates++;
1419 // }
1420 // }
1421 // }
1422 //
1423 // sphereBlockStart[sI] = nTotalCandidates;
1424 // sphereBlockSize[sI] = nCandidates;
1425 //
1426 // nTotalCandidates += nCandidates;
1427 // }
1428 //
1429 // // Needed for reverseDistribute
1436 //
1437 // // TODO: This is faster, but results in more vertices being referred
1438 // boolList sphereOverlapsCandidate(testCentres.size(), true);
1475 //
1481 //
1482 // labelListList sphereProcs(centres.size());
1483 //
1484 // // Working storage for assessing processors
1485 // DynamicList<label> tmpProcs;
1486 //
1487 // forAll(centres, sI)
1488 // {
1489 // tmpProcs.clear();
1490 //
1491 // // Extract the sub list of results for this point
1492 //
1493 // SubList<bool> sphereProcResults
1494 // (
1495 // sphereOverlapsCandidate,
1496 // sphereBlockSize[sI],
1497 // sphereBlockStart[sI]
1498 // );
1499 //
1500 // forAll(sphereProcResults, sPRI)
1501 // {
1502 // if (sphereProcResults[sPRI])
1503 // {
1504 // tmpProcs.append(toCandidateProc[sphereBlockStart[sI] + sPRI]);
1505 // }
1506 // }
1507 //
1508 // sphereProcs[sI] = tmpProcs;
1509 // }
1510 //
1511 // return sphereProcs;
1512 //}
1513 
1514 
1515 //Foam::labelListList Foam::backgroundMeshDecomposition::overlapProcessors
1516 //(
1517 // const point& cc,
1518 // const scalar rSqr
1519 //) const
1520 //{
1521 // DynamicList<label> toCandidateProc;
1522 // label sphereBlockStart(-1);
1523 // label sphereBlockSize(-1);
1524 //
1525 // label nCandidates = 0;
1526 //
1527 // forAll(allBackgroundMeshBounds_, proci)
1528 // {
1529 // // It is assumed that the sphere in question overlaps the source
1530 // // processor, so don't test it, unless includeOwnProcessor is true
1531 // if
1532 // (
1533 // (includeOwnProcessor || proci != Pstream::myProcNo())
1534 // && allBackgroundMeshBounds_[proci].overlaps(cc, rSqr)
1535 // )
1536 // {
1537 // toCandidateProc.append(proci);
1538 //
1539 // nCandidates++;
1540 // }
1541 // }
1542 //
1543 // sphereBlockSize = nCandidates;
1544 // nTotalCandidates += nCandidates;
1545 //
1546 // // Needed for reverseDistribute
1547 // label preDistributionToCandidateProcSize = toCandidateProc.size();
1548 //
1549 // autoPtr<mapDistribute> map(buildMap(toCandidateProc));
1550 //
1551 // map().distribute(testCentres);
1552 // map().distribute(testRadiusSqrs);
1553 //
1554 // // TODO: This is faster, but results in more vertices being referred
1556 // boolList sphereOverlapsCandidate(testCentres.size(), false);
1557 //
1558 // // Test candidate spheres on candidate processors
1559 // forAll(testCentres, sI)
1560 // {
1561 // const point& c = testCentres[sI];
1562 // const scalar rSqr = testRadiusSqrs[sI];
1563 //
1564 // const bool flagOverlap = bFTreePtr_().findNearest(c, rSqr).hit();
1565 //
1566 // if (flagOverlap)
1567 // {
1568 // //if (vertexOctree.findAnyOverlap(c, rSqr))
1573 //
1588 // sphereOverlapsCandidate[sI] = true;
1590 // }
1591 // }
1592 //
1593 // map().reverseDistribute
1594 // (
1595 // preDistributionToCandidateProcSize,
1596 // sphereOverlapsCandidate
1597 // );
1598 //
1599 // labelListList sphereProcs(centres.size());
1600 //
1601 // // Working storage for assessing processors
1602 // DynamicList<label> tmpProcs;
1603 //
1604 // forAll(centres, sI)
1605 // {
1606 // tmpProcs.clear();
1607 //
1608 // // Extract the sub list of results for this point
1609 //
1610 // SubList<bool> sphereProcResults
1611 // (
1612 // sphereOverlapsCandidate,
1613 // sphereBlockSize[sI],
1614 // sphereBlockStart[sI]
1615 // );
1616 //
1617 // forAll(sphereProcResults, sPRI)
1618 // {
1619 // if (sphereProcResults[sPRI])
1620 // {
1621 // tmpProcs.append(toCandidateProc[sphereBlockStart[sI] + sPRI]);
1622 // }
1623 // }
1624 //
1625 // sphereProcs[sI] = tmpProcs;
1626 // }
1627 //
1628 // return sphereProcs;
1629 //}
1630 
1631 
1632 // ************************************************************************* //
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.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
labelList processorNearestPosition(const List< point > &pts) const
What is the nearest processor to the given position?
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:319
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:417
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
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
labelList overlapProcessors(const point &centre, const scalar radiusSqr) const
bool positionOnThisProcessor(const point &pt) const
Is the given position inside the domain of this decomposition.
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
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:52
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
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:139
stressControl lookup("compactNormalStress") >> compactNormalStress
label nPoints
word timeName
Definition: getTimeIndex.H:3
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:91
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
PrimitivePatch< face, List, const pointField > bPatch
label readLabel(Istream &is)
Definition: label.H:64
dimensioned< scalar > magSqr(const dimensioned< Type > &)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
static const char nl
Definition: Ostream.H:262
defineTypeNameAndDebug(combustionModel, 0)
pointIndexHit findLine(const point &start, const point &end) const
Find nearest intersection of line between start and end, (exposing.
labelList f(nPoints)
bool overlapsThisProcessor(const treeBoundBox &box) const
Does the given box overlap the faces of the boundary of this.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
PrimitivePatch< face, SubList, const pointField & > primitivePatch
Addressing for a faceList slice.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
bool overlapsOtherProcessors(const point &centre, const scalar &radiusSqr) const
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:393
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:399
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
vector point
Point is a vector.
Definition: point.H:41
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
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,.
virtual Ostream & write(const token &)=0
Write next token to stream.
volScalarField & p
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.
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:1020
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end, (exposing.