trackingInverseDistanceCellCellStencil.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) 2017-2020 OpenCFD Ltd.
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify i
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
30 #include "OBJstream.H"
31 #include "Time.H"
32 #include "fvMeshSubset.H"
33 #include "globalIndex.H"
34 #include "oversetFvPatch.H"
36 #include "syncTools.H"
37 #include "treeBoundBoxList.H"
38 #include "voxelMeshSearch.H"
39 #include "dynamicOversetFvMesh.H"
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
43 namespace Foam
44 {
45 namespace cellCellStencils
46 {
47  defineTypeNameAndDebug(trackingInverseDistance, 0);
48  addToRunTimeSelectionTable(cellCellStencil, trackingInverseDistance, mesh);
49 }
50 }
51 
52 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
53 
55 (
56  const fvMesh& mesh,
57  const vector& smallVec,
58 
59  const boundBox& bb,
60  labelVector& nDivs,
62 
63  const labelList& cellMap,
64  labelList& patchCellTypes
65 )
66 {
67  // Mark all voxels that overlap the bounding box of any patch
68 
69  static bool hasWarned = false;
70 
71  const fvBoundaryMesh& pbm = mesh.boundary();
72 
73  patchTypes = patchCellType::OTHER;
74 
75  // Mark wall boundaries
76  forAll(pbm, patchi)
77  {
78  const fvPatch& fvp = pbm[patchi];
79  const labelList& fc = fvp.faceCells();
80 
81  if (!fvPatch::constraintType(fvp.type()))
82  {
83  //Info<< "Marking cells on proper patch " << fvp.name()
84  // << " with type " << fvp.type() << endl;
85  const polyPatch& pp = fvp.patch();
86  forAll(pp, i)
87  {
88  // Mark in overall patch types
89  patchCellTypes[cellMap[fc[i]]] = patchCellType::PATCH;
90 
91  // Mark in voxel mesh
92  boundBox faceBb(pp.points(), pp[i], false);
93  faceBb.min() -= smallVec;
94  faceBb.max() += smallVec;
95 
96  if (bb.overlaps(faceBb))
97  {
98  voxelMeshSearch::fill
99  (
100  patchTypes,
101  bb,
102  nDivs,
103  faceBb,
105  );
106  }
107  }
108  }
109  }
110 
111  // Override with overset boundaries
112  forAll(pbm, patchi)
113  {
114  const fvPatch& fvp = pbm[patchi];
115  const labelList& fc = fvp.faceCells();
116 
117  if (isA<oversetFvPatch>(fvp))
118  {
119  //Info<< "Marking cells on overset patch " << fvp.name() << endl;
120  const polyPatch& pp = fvp.patch();
121  const vectorField::subField faceCentres(pp.faceCentres());
122  forAll(pp, i)
123  {
124  // Mark in overall patch types
125  patchCellTypes[cellMap[fc[i]]] = patchCellType::OVERSET;
126 
127  // Mark in voxel mesh
128  boundBox faceBb(pp.points(), pp[i], false);
129  faceBb.min() -= smallVec;
130  faceBb.max() += smallVec;
131  if (!bb.contains(faceCentres[i]))
132  {
133  if (!hasWarned)
134  {
135  hasWarned = true;
136  WarningInFunction << "Patch " << pp.name()
137  << " : face at " << faceCentres[i]
138  << " is not inside search bounding box of"
139  << " voxel mesh " << bb << endl
140  << " Is your 'searchBox' specification correct?"
141  << endl;
142  }
143  }
144  else
145  {
146  // Test for having voxels already marked as patch
147  // -> voxel mesh is too coarse
148  if
149  (
150  voxelMeshSearch::overlaps
151  (
152  bb,
153  nDivs,
154  faceBb,
155  patchTypes,
156  static_cast<unsigned int>(patchCellType::PATCH)
157  )
158  )
159  {
160  // Determine number of voxels from number of cells
161  // in mesh
162  const labelVector& dim = mesh.geometricD();
163  forAll(dim, i)
164  {
165  if (dim[i] != -1)
166  {
167  nDivs[i] *= 2;
168  }
169  }
170 
171  Pout<< "Voxel mesh too coarse. Bounding box "
172  << faceBb
173  << " contains both non-overset and overset patches"
174  << ". Refining voxel mesh to " << nDivs << endl;
175 
176  return false;
177  }
178 
179  voxelMeshSearch::fill
180  (
181  patchTypes,
182  bb,
183  nDivs,
184  faceBb,
185  patchCellType::OVERSET
186  );
187  }
188  }
189  }
190  }
191 
192  return true;
193 }
194 
195 
197 (
198  PstreamBuffers& pBufs,
199  const List<treeBoundBoxList>& patchBb,
200  const List<labelVector>& patchDivisions,
201  const PtrList<PackedList<2>>& patchParts,
202 
203  const label srcI,
204  const label tgtI,
205  labelList& allCellTypes
206 ) const
207 {
208  const pointField& allPoints = mesh_.points();
209  const labelListList& allCellPoints = mesh_.cellPoints();
210 
211  const treeBoundBoxList& srcPatchBbs = patchBb[srcI];
212  const treeBoundBoxList& tgtPatchBbs = patchBb[tgtI];
213  const labelList& tgtCellMap = meshParts_[tgtI].cellMap();
214 
215  // 1. do processor-local src-tgt patch overlap
216  {
217  const treeBoundBox& srcPatchBb = srcPatchBbs[Pstream::myProcNo()];
218  const treeBoundBox& tgtPatchBb = tgtPatchBbs[Pstream::myProcNo()];
219 
220  if (srcPatchBb.overlaps(tgtPatchBb))
221  {
222  const PackedList<2>& srcPatchTypes = patchParts[srcI];
223  const labelVector& srcDivs = patchDivisions[srcI];
224 
225  forAll(tgtCellMap, tgtCelli)
226  {
227  label celli = tgtCellMap[tgtCelli];
228  boundBox cBb(allPoints, allCellPoints[celli], false);
229  cBb.min() -= smallVec_;
230  cBb.max() += smallVec_;
231 
232  if
233  (
234  voxelMeshSearch::overlaps
235  (
236  srcPatchBb,
237  srcDivs,
238  cBb,
239  srcPatchTypes,
240  static_cast<unsigned int>(patchCellType::PATCH)
241  )
242  )
243  {
244  allCellTypes[celli] = HOLE;
245  }
246  }
247  }
248  }
249 
250 
251  // 2. Send over srcMesh bits that overlap tgt and do calculation
252  pBufs.clear();
253  for (const int procI : Pstream::allProcs())
254  {
255  if (procI != Pstream::myProcNo())
256  {
257  const treeBoundBox& srcPatchBb = srcPatchBbs[Pstream::myProcNo()];
258  const treeBoundBox& tgtPatchBb = tgtPatchBbs[procI];
259 
260  if (srcPatchBb.overlaps(tgtPatchBb))
261  {
262  // Send over complete patch voxel map. Tbd: could
263  // subset
264  UOPstream os(procI, pBufs);
265  os << srcPatchBb << patchDivisions[srcI] << patchParts[srcI];
266  }
267  }
268  }
269  pBufs.finishedSends();
270  for (const int procI : Pstream::allProcs())
271  {
272  if (procI != Pstream::myProcNo())
273  {
274  const treeBoundBox& srcPatchBb = srcPatchBbs[procI];
275  const treeBoundBox& tgtPatchBb = tgtPatchBbs[Pstream::myProcNo()];
276 
277  if (srcPatchBb.overlaps(tgtPatchBb))
278  {
279  UIPstream is(procI, pBufs);
280  {
281  treeBoundBox receivedBb(is);
282  if (srcPatchBb != receivedBb)
283  {
285  << "proc:" << procI
286  << " srcPatchBb:" << srcPatchBb
287  << " receivedBb:" << receivedBb
288  << exit(FatalError);
289  }
290  }
291  const labelVector srcDivs(is);
292  const PackedList<2> srcPatchTypes(is);
293 
294  forAll(tgtCellMap, tgtCelli)
295  {
296  label celli = tgtCellMap[tgtCelli];
297  boundBox cBb(allPoints, allCellPoints[celli], false);
298  cBb.min() -= smallVec_;
299  cBb.max() += smallVec_;
300 
301  if
302  (
303  voxelMeshSearch::overlaps
304  (
305  srcPatchBb,
306  srcDivs,
307  cBb,
308  srcPatchTypes,
309  static_cast<unsigned int>(patchCellType::PATCH)
310  )
311  )
312  {
313  allCellTypes[celli] = HOLE;
314  }
315  }
316  }
317  }
318  }
319 }
320 
321 
323 (
324  PstreamBuffers& pBufs,
326  const PtrList<voxelMeshSearch>& meshSearches,
327  const labelList& allCellTypes,
328 
329  const label srcI,
330  const label tgtI,
331  labelListList& allStencil,
332  labelList& allDonor
333 ) const
334 {
335  const treeBoundBoxList& srcBbs = meshBb[srcI];
336  const treeBoundBoxList& tgtBbs = meshBb[tgtI];
337 
338  const fvMesh& srcMesh = meshParts_[srcI].subMesh();
339  const labelList& srcCellMap = meshParts_[srcI].cellMap();
340  const voxelMeshSearch& meshSearch = meshSearches[srcI];
341  const fvMesh& tgtMesh = meshParts_[tgtI].subMesh();
342  const pointField& tgtCc = tgtMesh.cellCentres();
343  const labelList& tgtCellMap = meshParts_[tgtI].cellMap();
344 
345  // 1. do processor-local src/tgt overlap
346  {
347  forAll(tgtCellMap, tgtCelli)
348  {
349  label srcCelli = meshSearch.findCell(tgtCc[tgtCelli]);
350  if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
351  {
352  label celli = tgtCellMap[tgtCelli];
353 
354  // TBD: check for multiple donors. Maybe better one? For
355  // now check 'nearer' mesh
356  if (betterDonor(tgtI, allDonor[celli], srcI))
357  {
358  allStencil[celli].setSize(1);
359  allStencil[celli][0] =
360  globalCells_.toGlobal(srcCellMap[srcCelli]);
361  allDonor[celli] = srcI;
362  }
363  }
364  }
365  }
366 
367 
368  // 2. Send over tgtMesh bits that overlap src and do calculation on
369  // srcMesh.
370 
371 
372  // (remote) processors where the tgt overlaps my src
373  DynamicList<label> tgtOverlapProcs(Pstream::nProcs());
374  // (remote) processors where the src overlaps my tgt
375  DynamicList<label> srcOverlapProcs(Pstream::nProcs());
376  for (const int procI : Pstream::allProcs())
377  {
378  if (procI != Pstream::myProcNo())
379  {
380  if (tgtBbs[procI].overlaps(srcBbs[Pstream::myProcNo()]))
381  {
382  tgtOverlapProcs.append(procI);
383  }
384  if (srcBbs[procI].overlaps(tgtBbs[Pstream::myProcNo()]))
385  {
386  srcOverlapProcs.append(procI);
387  }
388  }
389  }
390 
391 
392 
393  // Indices of tgtcells to send over to each processor
394  List<DynamicList<label>> tgtSendCells(Pstream::nProcs());
395  forAll(srcOverlapProcs, i)
396  {
397  label procI = srcOverlapProcs[i];
398  tgtSendCells[procI].reserve(tgtMesh.nCells()/srcOverlapProcs.size());
399  }
400 
401 
402  forAll(tgtCellMap, tgtCelli)
403  {
404  label celli = tgtCellMap[tgtCelli];
405  if (srcOverlapProcs.size())
406  {
407  treeBoundBox subBb(cellBb(mesh_, celli));
408  subBb.min() -= smallVec_;
409  subBb.max() += smallVec_;
410 
411  forAll(srcOverlapProcs, i)
412  {
413  label procI = srcOverlapProcs[i];
414  if (subBb.overlaps(srcBbs[procI]))
415  {
416  tgtSendCells[procI].append(tgtCelli);
417  }
418  }
419  }
420  }
421 
422  // Send target cell centres to overlapping processors
423  pBufs.clear();
424 
425  forAll(srcOverlapProcs, i)
426  {
427  label procI = srcOverlapProcs[i];
428  const labelList& cellIDs = tgtSendCells[procI];
429 
430  UOPstream os(procI, pBufs);
431  os << UIndirectList<point>(tgtCc, cellIDs);
432  }
433  pBufs.finishedSends();
434 
435  // Receive bits of target processors; find; send back
436  (void)srcMesh.tetBasePtIs();
437  forAll(tgtOverlapProcs, i)
438  {
439  label procI = tgtOverlapProcs[i];
440 
441  UIPstream is(procI, pBufs);
442  pointList samples(is);
443 
444  labelList donors(samples.size(), -1);
445  forAll(samples, sampleI)
446  {
447  label srcCelli = meshSearch.findCell(samples[sampleI]);
448  if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
449  {
450  donors[sampleI] = globalCells_.toGlobal(srcCellMap[srcCelli]);
451  }
452  }
453 
454  // Use same pStreamBuffers to send back.
455  UOPstream os(procI, pBufs);
456  os << donors;
457  }
458  pBufs.finishedSends();
459 
460  forAll(srcOverlapProcs, i)
461  {
462  label procI = srcOverlapProcs[i];
463  const labelList& cellIDs = tgtSendCells[procI];
464 
465  UIPstream is(procI, pBufs);
466  labelList donors(is);
467 
468  if (donors.size() != cellIDs.size())
469  {
470  FatalErrorInFunction<< "problem : cellIDs:" << cellIDs.size()
471  << " donors:" << donors.size() << abort(FatalError);
472  }
473 
474  forAll(donors, donorI)
475  {
476  label globalDonor = donors[donorI];
477 
478  if (globalDonor != -1)
479  {
480  label celli = tgtCellMap[cellIDs[donorI]];
481 
482  // TBD: check for multiple donors. Maybe better one?
483  if (betterDonor(tgtI, allDonor[celli], srcI))
484  {
485  allStencil[celli].setSize(1);
486  allStencil[celli][0] = globalDonor;
487  allDonor[celli] = srcI;
488  }
489  }
490  }
491  }
492 }
493 
494 
495 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
496 
497 Foam::cellCellStencils::trackingInverseDistance::trackingInverseDistance
498 (
499  const fvMesh& mesh,
500  const dictionary& dict,
501  const bool doUpdate
502 )
503 :
504  inverseDistance(mesh, dict, false),
505  globalCells_(mesh_.nCells())
506 {
507  // Initialise donor cell
508  globalDonor_.setSize(mesh_.nCells());
509  globalDonor_ = -1;
510 
511  // Initialise mesh partitions
512  const labelIOList& zoneID = this->zoneID();
513  label nZones = gMax(zoneID)+1;
514 
516  forAll(zoneID, celli)
517  {
518  nCellsPerZone[zoneID[celli]]++;
519  }
520  Pstream::listCombineGather(nCellsPerZone, plusEqOp<label>());
521  Pstream::listCombineScatter(nCellsPerZone);
522 
523  meshParts_.setSize(nZones);
524  forAll(meshParts_, zonei)
525  {
526  meshParts_.set
527  (
528  zonei,
529  new fvMeshSubset(mesh_, zonei, zoneID)
530  );
531 
532  // Trigger early evaluation of mesh dimension
533  // (in case there are locally zero cells in mesh)
534  (void)meshParts_[zonei].subMesh().nGeometricD();
535  }
536 
537 
538  // Print a bit
539  {
540  Info<< typeName << " : detected " << nZones
541  << " mesh regions" << endl;
542  Info<< incrIndent;
543  forAll(nCellsPerZone, zonei)
544  {
545  Info<< indent<< "zone:" << zonei
546  << " nCells:" << nCellsPerZone[zonei]
547  << endl;
548  }
549  Info<< decrIndent;
550  }
551 
552 
553  // Do geometry update
554  if (doUpdate)
555  {
556  update();
557  }
558 }
559 
560 
561 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
562 
564 {}
565 
566 
567 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
568 
570 {
571  DebugInfo<< FUNCTION_NAME << " : Start of analysis" << endl;
572 
573  scalar layerRelax(dict_.getOrDefault("layerRelax", 1.0));
574  const labelIOList& zoneID = this->zoneID();
575  label nZones = meshParts_.size();
576 
577  // Update stored mesh partitions for geometry changes
578  forAll(meshParts_, zonei)
579  {
580  pointField subPoints(mesh_.points(), meshParts_[zonei].pointMap());
581 
582  fvMesh& subMesh = meshParts_[zonei].subMesh();
583  subMesh.movePoints(subPoints);
584  }
585 
586  DebugInfo<< FUNCTION_NAME << " : Moved zone sub-meshes" << endl;
587 
588  // Calculate fast search structure for each zone
589  List<labelVector> searchBoxDivisions(dict_.lookup("searchBoxDivisions"));
590  if (searchBoxDivisions.size() != nZones)
591  {
592  FatalIOErrorInFunction(dict_) << "Number of searchBoxDivisions "
593  << searchBoxDivisions.size()
594  << " should equal the number of zones " << nZones
595  << exit(FatalIOError);
596  }
597 
598  const boundBox& allBb(mesh_.bounds());
599 
601 
602  // Determine zone meshes and bounding boxes
603  {
604  // Per processor, per zone the bounding box
606  procBb[Pstream::myProcNo()].setSize(nZones);
607 
608  forAll(meshParts_, zonei)
609  {
610  const fvMesh& subMesh = meshParts_[zonei].subMesh();
611 
612  if (subMesh.nPoints())
613  {
614  procBb[Pstream::myProcNo()][zonei] =
615  treeBoundBox(subMesh.points());
616  procBb[Pstream::myProcNo()][zonei].inflate(1e-6);
617  }
618  else
619  {
620  // No part of zone on this processor. Make up bb.
621  procBb[Pstream::myProcNo()][zonei] = treeBoundBox
622  (
623  allBb.min() - 2*allBb.span(),
624  allBb.min() - allBb.span()
625  );
626  procBb[Pstream::myProcNo()][zonei].inflate(1e-6);
627  }
628  }
629 
630  Pstream::gatherList(procBb);
631  Pstream::scatterList(procBb);
632 
633  // Move local bounding boxes to per-mesh indexing
634  forAll(meshBb, zonei)
635  {
636  treeBoundBoxList& bbs = meshBb[zonei];
637  bbs.setSize(Pstream::nProcs());
638  forAll(procBb, proci)
639  {
640  bbs[proci] = procBb[proci][zonei];
641  }
642  }
643  }
644 
645  DebugInfo<< FUNCTION_NAME << " : Calculated bounding boxes" << endl;
646 
647 
648  // Determine patch bounding boxes. These are either global and provided
649  // by the user or processor-local as a copy of the mesh bounding box.
650 
652  List<labelVector> patchDivisions(searchBoxDivisions);
653  PtrList<PackedList<2>> patchParts(nZones);
654  labelList allPatchTypes(mesh_.nCells(), OTHER);
655 
656  {
657  treeBoundBox globalPatchBb;
658  if (dict_.readIfPresent("searchBox", globalPatchBb))
659  {
660  // All processors, all zones have the same bounding box
661  patchBb = treeBoundBoxList(Pstream::nProcs(), globalPatchBb);
662  }
663  else
664  {
665  // Use the meshBb (differing per zone, per processor)
666  patchBb = meshBb;
667  }
668  }
669 
670  forAll(patchParts, zonei)
671  {
672  while (true)
673  {
674  patchParts.set
675  (
676  zonei,
677  new PackedList<2>(cmptProduct(patchDivisions[zonei]))
678  );
679 
680  bool ok = markBoundaries
681  (
682  meshParts_[zonei].subMesh(),
683  smallVec_,
684 
685  patchBb[zonei][Pstream::myProcNo()],
686  patchDivisions[zonei],
687  patchParts[zonei],
688 
689  meshParts_[zonei].cellMap(),
690  allPatchTypes
691  );
692 
693  //if (returnReduce(ok, andOp<bool>()))
694  if (ok)
695  {
696  break;
697  }
698  }
699  }
700  DebugInfo<< FUNCTION_NAME << " : Calculated boundary voxel meshes" << endl;
701 
702 
703  PtrList<voxelMeshSearch> meshSearches(meshParts_.size());
704  forAll(meshParts_, zonei)
705  {
706  meshSearches.set
707  (
708  zonei,
709  new voxelMeshSearch
710  (
711  meshParts_[zonei].subMesh(),
712  patchBb[zonei][Pstream::myProcNo()],
713  patchDivisions[zonei],
714  true
715  )
716  );
717  }
718 
719  DebugInfo<< FUNCTION_NAME << " : Constructed cell voxel meshes" << endl;
720 
721  PtrList<fvMesh> voxelMeshes;
722  if (debug)
723  {
724  voxelMeshes.setSize(meshSearches.size());
725  forAll(meshSearches, zonei)
726  {
727  IOobject meshIO
728  (
729  word("voxelMesh" + Foam::name(zonei)),
730  mesh_.time().timeName(),
731  mesh_,
733  );
734 
735  Pout<< "Constructing voxel mesh " << meshIO.objectPath() << endl;
736  voxelMeshes.set(zonei, meshSearches[zonei].makeMesh(meshIO));
737  }
738  }
739 
740 
741  if (debug)
742  {
743  forAll(patchParts, zonei)
744  {
745  const labelList voxelTypes(patchParts[zonei].values());
746  const fvMesh& vm = voxelMeshes[zonei];
748  (
749  createField<label>
750  (
751  vm,
752  "patchTypes",
753  voxelTypes
754  )
755  );
756  tfld().write();
757  }
758  }
759 
760  // Current best guess for cells
761  labelList allCellTypes(mesh_.nCells(), CALCULATED);
762  labelListList allStencil(mesh_.nCells());
763  // zoneID of donor
764  labelList allDonorID(mesh_.nCells(), -1);
765 
766  const globalIndex globalCells(mesh_.nCells());
767 
769 
770  DebugInfo<< FUNCTION_NAME << " : Allocated donor-cell structures" << endl;
771 
772  for (label srci = 0; srci < meshParts_.size()-1; srci++)
773  {
774  for (label tgti = srci+1; tgti < meshParts_.size(); tgti++)
775  {
776  markPatchesAsHoles
777  (
778  pBufs,
779 
780  patchBb,
781  patchDivisions,
782  patchParts,
783 
784  srci,
785  tgti,
786  allCellTypes
787  );
788  markPatchesAsHoles
789  (
790  pBufs,
791 
792  patchBb,
793  patchDivisions,
794  patchParts,
795 
796  tgti,
797  srci,
798  allCellTypes
799  );
800  }
801  }
802 
803  for (label srci = 0; srci < meshParts_.size()-1; srci++)
804  {
805  for (label tgti = srci+1; tgti < meshParts_.size(); tgti++)
806  {
807  markDonors
808  (
809  pBufs,
810  meshBb,
811  meshSearches,
812  allCellTypes, // to exclude hole donors
813 
814  tgti,
815  srci,
816  allStencil,
817  allDonorID
818  );
819  markDonors
820  (
821  pBufs,
822  meshBb,
823  meshSearches,
824  allCellTypes, // to exclude hole donors
825 
826  srci,
827  tgti,
828  allStencil,
829  allDonorID
830  );
831  }
832  }
833 
834  DebugInfo<< FUNCTION_NAME << " : Determined holes and donor-acceptors"
835  << endl;
836 
837  if (debug)
838  {
840  (
841  createField(mesh_, "allCellTypes", allCellTypes)
842  );
843  tfld().write();
844  }
845 
846 
847  // Use the patch types and weights to decide what to do
848  forAll(allPatchTypes, celli)
849  {
850  if (allCellTypes[celli] != HOLE)
851  {
852  switch (allPatchTypes[celli])
853  {
854  case OVERSET:
855  {
856  // Require interpolation. See if possible.
857 
858  if (allStencil[celli].size())
859  {
860  allCellTypes[celli] = INTERPOLATED;
861  }
862  else
863  {
864  allCellTypes[celli] = HOLE;
865  }
866  }
867  }
868  }
869  }
870  DebugInfo<< FUNCTION_NAME << " : Removed bad donors" << endl;
871 
872  if (debug)
873  {
875  (
876  createField(mesh_, "allCellTypes_patch", allCellTypes)
877  );
878  tfld().write();
879  }
880 
881 
882  // Check previous iteration cellTypes_ for any hole->calculated changes
883  // If so set the cell either to interpolated (if there are donors) or
884  // holes (if there are no donors). Note that any interpolated cell might
885  // still be overwritten by the flood filling
886  {
887  label nCalculated = 0;
888 
889  forAll(cellTypes_, celli)
890  {
891  if (allCellTypes[celli] == CALCULATED && cellTypes_[celli] == HOLE)
892  {
893  if (allStencil[celli].size() == 0)
894  {
895  // Reset to hole
896  allCellTypes[celli] = HOLE;
897  allStencil[celli].clear();
898  }
899  else
900  {
901  allCellTypes[celli] = INTERPOLATED;
902  nCalculated++;
903  }
904  }
905  }
906 
907  if (debug)
908  {
909  Pout<< "Detected " << nCalculated << " cells changing from hole"
910  << " to calculated. Changed to interpolated"
911  << endl;
912  }
913  }
914 
915 
916  // Mark unreachable bits
917  findHoles(globalCells_, mesh_, zoneID, allStencil, allCellTypes);
918  DebugInfo<< FUNCTION_NAME << " : Flood-filled holes" << endl;
919 
920  if (debug)
921  {
923  (
924  createField(mesh_, "allCellTypes_hole", allCellTypes)
925  );
926  tfld().write();
927  }
928 
929  // Add buffer interpolation layer(s) around holes
930  scalarField allWeight(mesh_.nCells(), Zero);
931  walkFront(layerRelax, allStencil, allCellTypes, allWeight);
932  DebugInfo<< FUNCTION_NAME << " : Implemented layer relaxation" << endl;
933 
934  if (debug)
935  {
937  (
938  createField(mesh_, "allCellTypes_front", allCellTypes)
939  );
940  tfld().write();
941  }
942 
943 
944  // Convert cell-cell addressing to stencil in compact notation
945 
946  cellTypes_.transfer(allCellTypes);
947  cellStencil_.setSize(mesh_.nCells());
948  cellInterpolationWeights_.setSize(mesh_.nCells());
949  DynamicList<label> interpolationCells;
950  forAll(cellTypes_, celli)
951  {
952  if (cellTypes_[celli] == INTERPOLATED)
953  {
954  cellStencil_[celli].transfer(allStencil[celli]);
955  cellInterpolationWeights_[celli].setSize(1);
956  cellInterpolationWeights_[celli][0] = 1.0;
957  interpolationCells.append(celli);
958  }
959  else
960  {
961  cellStencil_[celli].clear();
962  cellInterpolationWeights_[celli].clear();
963  }
964  }
965  interpolationCells_.transfer(interpolationCells);
966 
967  List<Map<label>> compactMap;
968  cellInterpolationMap_.reset
969  (
970  new mapDistribute
971  (
972  globalCells,
973  cellStencil_,
974  compactMap
975  )
976  );
977  cellInterpolationWeight_.transfer(allWeight);
978  //cellInterpolationWeight_.correctBoundaryConditions();
980  <
983  >(cellInterpolationWeight_.boundaryFieldRef(), false);
984 
985 
986  if (debug & 2)
987  {
988  // Dump stencil
989  mkDir(mesh_.time().timePath());
990  OBJstream str(mesh_.time().timePath()/"injectionStencil.obj");
991  Pout<< typeName << " : dumping injectionStencil to "
992  << str.name() << endl;
993  pointField cc(mesh_.cellCentres());
994  cellInterpolationMap().distribute(cc);
995 
996  forAll(cellStencil_, celli)
997  {
998  const labelList& slots = cellStencil_[celli];
999  if (slots.size())
1000  {
1001  const point& accCc = mesh_.cellCentres()[celli];
1002  forAll(slots, i)
1003  {
1004  const point& donorCc = cc[slots[i]];
1005  str.write(linePointRef(accCc, 0.1*accCc+0.9*donorCc));
1006  }
1007  }
1008  }
1009  }
1010 
1011  DebugInfo<< FUNCTION_NAME << " : Transferred donor to stencil" << endl;
1012 
1013 
1014  // Extend stencil to get inverse distance weighted neighbours
1015  createStencil(globalCells);
1016  DebugInfo<< FUNCTION_NAME << " : Extended stencil" << endl;
1017 
1018 
1019  if (debug & 2)
1020  {
1021  // Dump weight
1022  cellInterpolationWeight_.instance() = mesh_.time().timeName();
1023  cellInterpolationWeight_.write();
1024 
1025  // Dump cell types
1026  volScalarField volTypes
1027  (
1028  IOobject
1029  (
1030  "cellTypes",
1031  mesh_.time().timeName(),
1032  mesh_,
1035  false
1036  ),
1037  mesh_,
1039  zeroGradientFvPatchScalarField::typeName
1040  );
1041 
1042  forAll(volTypes.internalField(), celli)
1043  {
1044  volTypes[celli] = cellTypes_[celli];
1045  }
1046  //volTypes.correctBoundaryConditions();
1048  <
1051  >(volTypes.boundaryFieldRef(), false);
1052  volTypes.write();
1053 
1054  // Dump stencil
1055  mkDir(mesh_.time().timePath());
1056  OBJstream str(mesh_.time().timePath()/"stencil.obj");
1057  Pout<< typeName << " : dumping to " << str.name() << endl;
1058  pointField cc(mesh_.cellCentres());
1059  cellInterpolationMap().distribute(cc);
1060 
1061  forAll(cellStencil_, celli)
1062  {
1063  const labelList& slots = cellStencil_[celli];
1064  if (slots.size())
1065  {
1066  const point& accCc = mesh_.cellCentres()[celli];
1067  forAll(slots, i)
1068  {
1069  const point& donorCc = cc[slots[i]];
1070  str.write(linePointRef(accCc, 0.1*accCc+0.9*donorCc));
1071  }
1072  }
1073  }
1074  }
1075 
1076 
1077  // Print some stats
1078  {
1079  labelList nCells(count(3, cellTypes_));
1080 
1081  label nLocal = 0;
1082  label nMixed = 0;
1083  label nRemote = 0;
1084  forAll(interpolationCells_, i)
1085  {
1086  label celli = interpolationCells_[i];
1087  const labelList& slots = cellStencil_[celli];
1088 
1089  bool hasLocal = false;
1090  bool hasRemote = false;
1091 
1092  forAll(slots, sloti)
1093  {
1094  if (slots[sloti] >= mesh_.nCells())
1095  {
1096  hasRemote = true;
1097  }
1098  else
1099  {
1100  hasLocal = true;
1101  }
1102  }
1103 
1104  if (hasRemote)
1105  {
1106  if (!hasLocal)
1107  {
1108  nRemote++;
1109  }
1110  else
1111  {
1112  nMixed++;
1113  }
1114  }
1115  else if (hasLocal)
1116  {
1117  nLocal++;
1118  }
1119  }
1120  reduce(nLocal, sumOp<label>());
1121  reduce(nMixed, sumOp<label>());
1122  reduce(nRemote, sumOp<label>());
1123 
1124  Info<< "Overset analysis : nCells : "
1125  << returnReduce(cellTypes_.size(), sumOp<label>()) << nl
1126  << incrIndent
1127  << indent << "calculated : " << nCells[CALCULATED] << nl
1128  << indent << "interpolated : " << nCells[INTERPOLATED]
1129  << " (interpolated from local:" << nLocal
1130  << " mixed local/remote:" << nMixed
1131  << " remote:" << nRemote << ")" << nl
1132  << indent << "hole : " << nCells[HOLE] << nl
1133  << decrIndent << endl;
1134  }
1135  DebugInfo<< FUNCTION_NAME << " : Finished analysis" << endl;
1136 
1137  // Tbd: detect if anything changed. Most likely it did!
1138  return true;
1139 }
1140 
1141 
1142 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::PrimitivePatch::points
const Field< point_type > & points() const noexcept
Return reference to global points.
Definition: PrimitivePatch.H:299
nZones
label nZones
Definition: interpolatedFaces.H:24
Foam::IOobject::NO_WRITE
Definition: IOobject.H:195
Foam::polyMesh::points
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1069
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
Foam::treeBoundBox::overlaps
bool overlaps(const boundBox &bb) const
Overlaps other bounding box?
Definition: boundBoxI.H:221
Foam::dynamicOversetFvMesh::correctBoundaryConditions
static void correctBoundaryConditions(typename GeoField::Boundary &bfld, const bool typeOnly)
Correct boundary conditions of certain type (typeOnly = true)
Definition: dynamicOversetFvMeshTemplates.C:117
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::meshSearch
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search.
Definition: meshSearch.H:60
Foam::UOPstream
Output inter-processor communications stream operating on external buffer.
Definition: UOPstream.H:57
Foam::cellCellStencils::inverseDistance
Inverse-distance-weighted interpolation stencil.
Definition: inverseDistanceCellCellStencil.H:67
Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
void markPatchesAsHoles(PstreamBuffers &pBufs, const List< treeBoundBoxList > &patchBb, const List< labelVector > &patchDivisions, const PtrList< PackedList< 2 >> &patchParts, const label srcI, const label tgtI, labelList &allCellTypes) const
Mark all cells overlapping (a voxel covered by) a src patch.
Definition: trackingInverseDistanceCellCellStencil.C:197
Foam::returnReduce
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Definition: PstreamReduceOps.H:94
Foam::OBJstream
OFstream that keeps track of vertices.
Definition: OBJstream.H:58
update
mesh update()
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::fvMeshSubset
Given the original mesh and the list of selected cells, it creates the mesh consisting only of the de...
Definition: fvMeshSubset.H:73
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::DynamicList< label >
Foam::HashTableOps::values
List< T > values(const HashTable< T, Key, Hash > &tbl, const bool doSort=false)
List of values from HashTable, optionally sorted.
Definition: HashOps.H:149
Foam::cellCellStencils::trackingInverseDistance::~trackingInverseDistance
virtual ~trackingInverseDistance()
Destructor.
Definition: trackingInverseDistanceCellCellStencil.C:563
fvMeshSubset.H
Foam::Pstream::scatterList
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
Definition: gatherScatterList.C:215
globalIndex.H
Foam::treeBoundBox
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:86
Foam::PstreamBuffers
Buffers for inter-processor communications streams (UOPstream, UIPstream).
Definition: PstreamBuffers.H:88
voxelMeshSearch.H
Foam::cmptProduct
Cmpt cmptProduct(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:596
dynamicOversetFvMesh.H
Foam::List::append
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
Foam::FatalIOError
IOerror FatalIOError
Foam::fvMesh::movePoints
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: fvMesh.C:865
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::PtrList::set
const T * set(const label i) const
Return const pointer to element (can be nullptr),.
Definition: PtrList.H:138
Foam::boundBox::max
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
meshBb
List< treeBoundBox > meshBb(1, treeBoundBox(boundBox(coarseMesh.points(), false)).extend(rndGen, 1e-3))
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::PstreamBuffers::clear
void clear()
Reset (clear) individual buffers and reset state.
Definition: PstreamBuffers.C:125
Foam::meshSearch::findCell
label findCell(const point &location, const label seedCelli=-1, const bool useTreeSearch=true) const
Find cell containing location.
Definition: meshSearch.C:780
Foam::cellCellStencils::trackingInverseDistance::markDonors
void markDonors(PstreamBuffers &pBufs, const List< treeBoundBoxList > &meshBb, const PtrList< voxelMeshSearch > &meshSearches, const labelList &allCellTypes, const label srcI, const label tgtI, labelListList &allStencil, labelList &allDonor) const
Definition: trackingInverseDistanceCellCellStencil.C:323
Foam::incrIndent
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:346
syncTools.H
Foam::cellCellStencils::trackingInverseDistance::update
virtual bool update()
Update stencils. Return false if nothing changed.
Definition: trackingInverseDistanceCellCellStencil.C:569
Foam::boundBox::min
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:91
Foam::sumOp
Definition: ops.H:213
Foam::primitiveMesh::nPoints
label nPoints() const noexcept
Number of mesh points.
Definition: primitiveMeshI.H:37
Foam::cellCellStencils::trackingInverseDistance::markBoundaries
static bool markBoundaries(const fvMesh &mesh, const vector &smallVec, const boundBox &bb, labelVector &nDivs, PackedList< 2 > &patchTypes, const labelList &cellMap, labelList &patchCellTypes)
Mark voxels of patchTypes with type of patch face.
Definition: trackingInverseDistanceCellCellStencil.C:55
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
patchTypes
wordList patchTypes(nPatches)
Foam::fvBoundaryMesh
Foam::fvBoundaryMesh.
Definition: fvBoundaryMesh.H:57
Foam::SubField
SubField is a Field obtained as a section of another Field, without its own allocation....
Definition: Field.H:64
Foam::oversetFvPatchField::write
virtual void write(Ostream &os) const
Write.
Definition: oversetFvPatchField.C:222
Foam::boundBox::overlaps
bool overlaps(const boundBox &bb) const
Overlaps/touches boundingBox?
Definition: boundBoxI.H:221
trackingInverseDistanceCellCellStencil.H
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:51
Foam::primitiveMesh::nCells
label nCells() const noexcept
Number of mesh cells.
Definition: primitiveMeshI.H:96
Foam::Field< vector >
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::polyPatch
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:68
Foam::OSstream::name
virtual const fileName & name() const
Get the name of the stream.
Definition: OSstream.H:107
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:222
Foam::DynamicList::append
DynamicList< T, SizeMin > & append(const T &val)
Append an element to the end of this list.
Definition: DynamicListI.H:511
Foam::mapDistribute
Class containing processor-to-processor mapping information.
Definition: mapDistribute.H:163
Foam::List::transfer
void transfer(List< T > &list)
Definition: List.C:456
Foam::fvPatch
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:65
Foam::PtrList::setSize
void setSize(const label newLen)
Same as resize()
Definition: PtrList.H:151
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:42
samples
scalarField samples(nIntervals, Zero)
Foam::PtrList
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: List.H:59
Foam::volScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
Foam::PstreamBuffers::finishedSends
void finishedSends(const bool block=true)
Mark all sends as having been done.
Definition: PstreamBuffers.C:73
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
os
OBJstream os(runTime.globalPath()/outputName)
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
zoneID
const labelIOList & zoneID
Definition: interpolatedFaces.H:22
Foam::Ostream::write
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:85
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::decrIndent
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:353
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::globalIndex
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
Foam::indent
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:339
nCellsPerZone
labelList nCellsPerZone(nZones, Zero)
treeBoundBoxList.H
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::cellCellStencils::addToRunTimeSelectionTable
addToRunTimeSelectionTable(cellCellStencil, cellVolumeWeight, mesh)
Foam::boundBox::contains
bool contains(const point &pt) const
Contains point? (inside or on edge)
Definition: boundBoxI.H:271
Time.H
Foam::primitiveMesh::cellCentres
const vectorField & cellCentres() const
Definition: primitiveMeshCellCentresAndVols.C:84
Foam::linePointRef
line< point, const point & > linePointRef
A line using referred points.
Definition: linePointRef.H:47
Foam::voxelMeshSearch
Fast, non-parallel searching in mesh without use of octree.
Definition: voxelMeshSearch.H:57
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::Pstream::gatherList
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
Definition: gatherScatterList.C:52
Foam::UPstream::commsTypes::nonBlocking
DebugInfo
#define DebugInfo
Report an information message using Foam::Info.
Definition: messageStream.H:382
Foam::UPstream::myProcNo
static int myProcNo(const label communicator=worldComm)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:463
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::polyPatch::faceCentres
const vectorField::subField faceCentres() const
Return face centres.
Definition: polyPatch.C:321
Foam::VectorSpace::max
static const Form max
Definition: VectorSpace.H:117
Foam::BitOps::count
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
Definition: BitOps.H:77
Foam::Vector< scalar >
Foam::List< label >
Foam::PackedList< 2 >
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::IOList< label >
Foam::List::clear
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
oversetFvPatch.H
FUNCTION_NAME
#define FUNCTION_NAME
Definition: messageStream.H:295
Foam::oversetFvPatchField
Boundary condition for use on overset patches. To be run in combination with special dynamicFvMesh ty...
Definition: oversetFvPatchField.H:56
Foam::fvPatch::faceCells
virtual const labelUList & faceCells() const
Return faceCells.
Definition: fvPatch.C:113
Foam::boundBox
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:63
Foam::fvPatch::patch
const polyPatch & patch() const
Return the polyPatch.
Definition: fvPatch.H:161
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
FatalIOErrorInFunction
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:473
Foam::DelaunayMeshTools::allPoints
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
Foam::patchIdentifier::name
const word & name() const noexcept
The patch name.
Definition: patchIdentifier.H:135
Foam::cellCellStencils::defineTypeNameAndDebug
defineTypeNameAndDebug(cellVolumeWeight, 0)
Foam::plusEqOp
Definition: ops.H:72
Foam::PATCH
PDRpatchDef PATCH
Definition: PDRpatchDef.H:112
Foam::UIPstream
Input inter-processor communications stream operating on external buffer.
Definition: UIPstream.H:56
Foam::GeometricField< scalar, fvPatchField, volMesh >
Foam::IOobject::NO_READ
Definition: IOobject.H:188
Foam::treeBoundBoxList
List< treeBoundBox > treeBoundBoxList
List of bounding boxes.
Definition: treeBoundBoxList.H:46
Foam::mkDir
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
Definition: MSwindows.C:507
zeroGradientFvPatchFields.H
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::gMax
Type gMax(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:592
OBJstream.H
Foam::UPstream::nProcs
static label nProcs(const label communicator=worldComm)
Number of processes in parallel run, and 1 for serial run.
Definition: UPstream.H:445
Foam::dimless
const dimensionSet dimless
Dimensionless.
Definition: dimensionSets.C:189
Foam::polyMesh::tetBasePtIs
const labelIOList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:892