InteractionLists.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-2022 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 
26 #include "InteractionLists.H"
28 #include "indexedOctree.H"
29 #include "treeDataFace.H"
30 #include "treeDataCell.H"
31 #include "volFields.H"
32 #include "meshTools.H"
33 
34 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
35 
36 template<class ParticleType>
38 {
39  Info<< "Building InteractionLists with interaction distance "
40  << maxDistance_ << endl;
41 
42  const vector interactionVec = maxDistance_*vector::one;
43 
44  treeBoundBox procBb(treeBoundBox(mesh_.points()));
45 
46  treeBoundBox extendedProcBb
47  (
48  procBb.min() - interactionVec,
49  procBb.max() + interactionVec
50  );
51 
52  treeBoundBoxList allExtendedProcBbs(Pstream::nProcs());
53 
54  allExtendedProcBbs[Pstream::myProcNo()] = extendedProcBb;
55 
56  Pstream::gatherList(allExtendedProcBbs);
57 
58  Pstream::scatterList(allExtendedProcBbs);
59 
60  List<treeBoundBox> extendedProcBbsInRange;
61  List<label> extendedProcBbsTransformIndex;
62  List<label> extendedProcBbsOrigProc;
63 
64  findExtendedProcBbsInRange
65  (
66  procBb,
67  allExtendedProcBbs,
68  mesh_.globalData().globalTransforms(),
69  extendedProcBbsInRange,
70  extendedProcBbsTransformIndex,
71  extendedProcBbsOrigProc
72  );
73 
74  treeBoundBoxList cellBbs(mesh_.nCells());
75 
76  forAll(cellBbs, celli)
77  {
78  cellBbs[celli] = treeBoundBox
79  (
80  mesh_.cells()[celli].points
81  (
82  mesh_.faces(),
83  mesh_.points()
84  )
85  );
86  }
87 
88  const globalIndexAndTransform& globalTransforms =
89  mesh_.globalData().globalTransforms();
90 
91  // Recording which cells are in range of an extended boundBox, as
92  // only these cells will need to be tested to determine which
93  // referred cells that they interact with.
94  PackedBoolList cellInRangeOfCoupledPatch(mesh_.nCells(), false);
95 
96  // IAndT: index (=local cell index) and transform (from
97  // globalIndexAndTransform)
98  DynamicList<labelPair> cellIAndTToExchange;
99 
100  DynamicList<treeBoundBox> cellBbsToExchange;
101 
102  DynamicList<label> procToDistributeCellTo;
103 
104  forAll(extendedProcBbsInRange, ePBIRI)
105  {
106  const treeBoundBox& otherExtendedProcBb =
107  extendedProcBbsInRange[ePBIRI];
108 
109  label transformIndex = extendedProcBbsTransformIndex[ePBIRI];
110 
111  label origProc = extendedProcBbsOrigProc[ePBIRI];
112 
113  forAll(cellBbs, celli)
114  {
115  const treeBoundBox& cellBb = cellBbs[celli];
116 
117  if (cellBb.overlaps(otherExtendedProcBb))
118  {
119  // This cell is in range of the Bb of the other
120  // processor Bb, and so needs to be referred to it
121 
122  cellInRangeOfCoupledPatch[celli] = true;
123 
124  cellIAndTToExchange.append
125  (
126  globalTransforms.encode(celli, transformIndex)
127  );
128 
129  cellBbsToExchange.append(cellBb);
130 
131  procToDistributeCellTo.append(origProc);
132  }
133  }
134  }
135 
136  buildMap(cellMapPtr_, procToDistributeCellTo);
137 
138  // Needed for reverseDistribute
139  label preDistributionCellMapSize = procToDistributeCellTo.size();
140 
141  cellMap().distribute(cellBbsToExchange);
142 
143  cellMap().distribute(cellIAndTToExchange);
144 
145  // Determine labelList specifying only cells that are in range of
146  // a coupled boundary to build an octree limited to these cells.
147  DynamicList<label> coupledPatchRangeCells;
148 
149  forAll(cellInRangeOfCoupledPatch, celli)
150  {
151  if (cellInRangeOfCoupledPatch[celli])
152  {
153  coupledPatchRangeCells.append(celli);
154  }
155  }
156 
157  treeBoundBox procBbRndExt
158  (
159  treeBoundBox(mesh_.points()).extend(1e-4)
160  );
161 
162  indexedOctree<treeDataCell> coupledPatchRangeTree
163  (
164  treeDataCell
165  (
166  true, // Cache cell bb
167  mesh_,
168  coupledPatchRangeCells, // Subset of mesh
169  polyMesh::CELL_TETS // Consistent with tracking
170  ),
171  procBbRndExt,
172  8, // maxLevel,
173  10, // leafSize,
174  100.0 // duplicity
175  );
176 
177  ril_.setSize(cellBbsToExchange.size());
178 
179  // This needs to be a boolList, not PackedBoolList if
180  // reverseDistribute is called.
181  boolList cellBbRequiredByAnyCell(cellBbsToExchange.size(), false);
182 
183  Info<< " Building referred interaction lists" << endl;
184 
185  forAll(cellBbsToExchange, bbI)
186  {
187  const labelPair& ciat = cellIAndTToExchange[bbI];
188 
189  const transformer& transform = globalTransforms.transform
190  (
191  globalTransforms.transformIndex(ciat)
192  );
193 
194  treeBoundBox tempTransformedBb
195  (
196  transform.invTransformPosition(cellBbsToExchange[bbI].points())()
197  );
198 
199  treeBoundBox extendedBb
200  (
201  tempTransformedBb.min() - interactionVec,
202  tempTransformedBb.max() + interactionVec
203  );
204 
205  // Find all elements intersecting box.
206  labelList interactingElems
207  (
208  coupledPatchRangeTree.findBox(extendedBb)
209  );
210 
211  if (!interactingElems.empty())
212  {
213  cellBbRequiredByAnyCell[bbI] = true;
214  }
215 
216  ril_[bbI].setSize(interactingElems.size(), -1);
217 
218  forAll(interactingElems, i)
219  {
220  label elemI = interactingElems[i];
221 
222  // Here, a more detailed geometric test could be applied,
223  // i.e. a more accurate bounding volume like a OBB or
224  // convex hull or an exact geometrical test.
225 
226  label c = coupledPatchRangeTree.shapes().cellLabels()[elemI];
227 
228  ril_[bbI][i] = c;
229  }
230  }
231 
232  // Perform subset of ril_, to remove any referred cells that do
233  // not interact. They will not be sent from originating
234  // processors. This assumes that the disappearance of the cell
235  // from the sending list of the source processor, simply removes
236  // the referred cell from the ril_, all of the subsequent indices
237  // shuffle down one, but the order and structure is preserved,
238  // i.e. it, is as if the cell had never been referred in the first
239  // place.
240 
241  inplaceSubset(cellBbRequiredByAnyCell, ril_);
242 
243  // Send information about which cells are actually required back
244  // to originating processors.
245 
246  // At this point, cellBbsToExchange does not need to be maintained
247  // or distributed as it is not longer needed.
248 
249  cellBbsToExchange.setSize(0);
250 
251  cellMap().reverseDistribute
252  (
253  preDistributionCellMapSize,
254  cellBbRequiredByAnyCell
255  );
256 
257  cellMap().reverseDistribute
258  (
259  preDistributionCellMapSize,
260  cellIAndTToExchange
261  );
262 
263  // Perform ordering of cells to send, this invalidates the
264  // previous value of preDistributionCellMapSize.
265 
266  preDistributionCellMapSize = -1;
267 
268  // Move all of the used cells to refer to the start of the list
269  // and truncate it
270 
271  inplaceSubset(cellBbRequiredByAnyCell, cellIAndTToExchange);
272 
273  inplaceSubset(cellBbRequiredByAnyCell, procToDistributeCellTo);
274 
275  preDistributionCellMapSize = procToDistributeCellTo.size();
276 
277  // Rebuild distributionMap with only required referred cells
278  buildMap(cellMapPtr_, procToDistributeCellTo);
279 
280  // Store cellIndexAndTransformToDistribute
281  cellIndexAndTransformToDistribute_.transfer(cellIAndTToExchange);
282 
283  // Determine inverse addressing for referred cells
284 
285  rilInverse_.setSize(mesh_.nCells());
286 
287  // Temporary Dynamic lists for accumulation
288  List<DynamicList<label>> rilInverseTemp(rilInverse_.size());
289 
290  // Loop over all referred cells
291  forAll(ril_, refCelli)
292  {
293  const labelList& realCells = ril_[refCelli];
294 
295  // Loop over all real cells in that the referred cell is to
296  // supply interactions to and record the index of this
297  // referred cell in the real cells entry in rilInverse
298 
299  forAll(realCells, realCelli)
300  {
301  rilInverseTemp[realCells[realCelli]].append(refCelli);
302  }
303  }
304 
305  forAll(rilInverse_, celli)
306  {
307  rilInverse_[celli].transfer(rilInverseTemp[celli]);
308  }
309 
310  // Determine which wall faces to refer
311 
312  // The referring of wall patch data relies on patches having the same
313  // index on each processor.
314  mesh_.boundaryMesh().checkParallelSync(true);
315 
316  // Determine the index of all of the wall faces on this processor
317  DynamicList<label> localWallFaces;
318 
319  forAll(mesh_.boundaryMesh(), patchi)
320  {
321  const polyPatch& patch = mesh_.boundaryMesh()[patchi];
322 
323  if (isA<wallPolyPatch>(patch))
324  {
325  localWallFaces.append(identity(patch.size()) + patch.start());
326  }
327  }
328 
329  treeBoundBoxList wallFaceBbs(localWallFaces.size());
330 
331  forAll(wallFaceBbs, i)
332  {
333  wallFaceBbs[i] = treeBoundBox
334  (
335  mesh_.faces()[localWallFaces[i]].points(mesh_.points())
336  );
337  }
338 
339  // IAndT: index and transform
340  DynamicList<labelPair> wallFaceIAndTToExchange;
341 
342  DynamicList<treeBoundBox> wallFaceBbsToExchange;
343 
344  DynamicList<label> procToDistributeWallFaceTo;
345 
346  forAll(extendedProcBbsInRange, ePBIRI)
347  {
348  const treeBoundBox& otherExtendedProcBb =
349  extendedProcBbsInRange[ePBIRI];
350 
351  label transformIndex = extendedProcBbsTransformIndex[ePBIRI];
352 
353  label origProc = extendedProcBbsOrigProc[ePBIRI];
354 
355  forAll(wallFaceBbs, i)
356  {
357  const treeBoundBox& wallFaceBb = wallFaceBbs[i];
358 
359  if (wallFaceBb.overlaps(otherExtendedProcBb))
360  {
361  // This wall face is in range of the Bb of the other
362  // processor Bb, and so needs to be referred to it
363 
364  label wallFacei = localWallFaces[i];
365 
366  wallFaceIAndTToExchange.append
367  (
368  globalTransforms.encode(wallFacei, transformIndex)
369  );
370 
371  wallFaceBbsToExchange.append(wallFaceBb);
372 
373  procToDistributeWallFaceTo.append(origProc);
374  }
375  }
376  }
377 
378  buildMap(wallFaceMapPtr_, procToDistributeWallFaceTo);
379 
380  // Needed for reverseDistribute
381  label preDistributionWallFaceMapSize = procToDistributeWallFaceTo.size();
382 
383  wallFaceMap().distribute(wallFaceBbsToExchange);
384 
385  wallFaceMap().distribute(wallFaceIAndTToExchange);
386 
387  indexedOctree<treeDataCell> allCellsTree
388  (
389  treeDataCell(true, mesh_, polyMesh::CELL_TETS),
390  procBbRndExt,
391  8, // maxLevel,
392  10, // leafSize,
393  100.0 // duplicity
394  );
395 
396  rwfil_.setSize(wallFaceBbsToExchange.size());
397 
398  // This needs to be a boolList, not PackedBoolList if
399  // reverseDistribute is called.
400  boolList wallFaceBbRequiredByAnyCell(wallFaceBbsToExchange.size(), false);
401 
402  forAll(wallFaceBbsToExchange, bbI)
403  {
404  const labelPair& wfiat = wallFaceIAndTToExchange[bbI];
405 
406  const transformer& transform = globalTransforms.transform
407  (
408  globalTransforms.transformIndex(wfiat)
409  );
410 
411  treeBoundBox tempTransformedBb
412  (
413  transform.invTransformPosition
414  (
415  wallFaceBbsToExchange[bbI].points()
416  )()
417  );
418 
419  treeBoundBox extendedBb
420  (
421  tempTransformedBb.min() - interactionVec,
422  tempTransformedBb.max() + interactionVec
423  );
424 
425  // Find all elements intersecting box.
426  labelList interactingElems
427  (
428  coupledPatchRangeTree.findBox(extendedBb)
429  );
430 
431  if (!interactingElems.empty())
432  {
433  wallFaceBbRequiredByAnyCell[bbI] = true;
434  }
435 
436  rwfil_[bbI].setSize(interactingElems.size(), -1);
437 
438  forAll(interactingElems, i)
439  {
440  label elemI = interactingElems[i];
441 
442  // Here, a more detailed geometric test could be applied,
443  // i.e. a more accurate bounding volume like a OBB or
444  // convex hull or an exact geometrical test.
445 
446  label c = coupledPatchRangeTree.shapes().cellLabels()[elemI];
447 
448  rwfil_[bbI][i] = c;
449  }
450  }
451 
452  // Perform subset of rwfil_, to remove any referred wallFaces that do
453  // not interact. They will not be sent from originating
454  // processors. This assumes that the disappearance of the wallFace
455  // from the sending list of the source processor, simply removes
456  // the referred wallFace from the rwfil_, all of the subsequent indices
457  // shuffle down one, but the order and structure is preserved,
458  // i.e. it, is as if the wallFace had never been referred in the first
459  // place.
460 
461  inplaceSubset(wallFaceBbRequiredByAnyCell, rwfil_);
462 
463  // Send information about which wallFaces are actually required
464  // back to originating processors.
465 
466  // At this point, wallFaceBbsToExchange does not need to be
467  // maintained or distributed as it is not longer needed.
468 
469  wallFaceBbsToExchange.setSize(0);
470 
471  wallFaceMap().reverseDistribute
472  (
473  preDistributionWallFaceMapSize,
474  wallFaceBbRequiredByAnyCell
475  );
476 
477  wallFaceMap().reverseDistribute
478  (
479  preDistributionWallFaceMapSize,
480  wallFaceIAndTToExchange
481  );
482 
483  // Perform ordering of wallFaces to send, this invalidates the
484  // previous value of preDistributionWallFaceMapSize.
485 
486  preDistributionWallFaceMapSize = -1;
487 
488  // Move all of the used wallFaces to refer to the start of the
489  // list and truncate it
490 
491  inplaceSubset(wallFaceBbRequiredByAnyCell, wallFaceIAndTToExchange);
492 
493  inplaceSubset(wallFaceBbRequiredByAnyCell, procToDistributeWallFaceTo);
494 
495  preDistributionWallFaceMapSize = procToDistributeWallFaceTo.size();
496 
497  // Rebuild distributionMap with only required referred wallFaces
498  buildMap(wallFaceMapPtr_, procToDistributeWallFaceTo);
499 
500  // Store wallFaceIndexAndTransformToDistribute
501  wallFaceIndexAndTransformToDistribute_.transfer(wallFaceIAndTToExchange);
502 
503  // Determine inverse addressing for referred wallFaces
504 
505  rwfilInverse_.setSize(mesh_.nCells());
506 
507  // Temporary Dynamic lists for accumulation
508  List<DynamicList<label>> rwfilInverseTemp(rwfilInverse_.size());
509 
510  // Loop over all referred wall faces
511  forAll(rwfil_, refWallFacei)
512  {
513  const labelList& realCells = rwfil_[refWallFacei];
514 
515  // Loop over all real cells in that the referred wall face is
516  // to supply interactions to and record the index of this
517  // referred wall face in the real cells entry in rwfilInverse
518 
519  forAll(realCells, realCelli)
520  {
521  rwfilInverseTemp[realCells[realCelli]].append(refWallFacei);
522  }
523  }
524 
525  forAll(rwfilInverse_, celli)
526  {
527  rwfilInverse_[celli].transfer(rwfilInverseTemp[celli]);
528  }
529 
530  // Refer wall faces to the appropriate processor
531  referredWallFaces_.setSize(wallFaceIndexAndTransformToDistribute_.size());
532 
533  forAll(referredWallFaces_, rWFI)
534  {
535  const labelPair& wfiat = wallFaceIndexAndTransformToDistribute_[rWFI];
536 
537  label wallFaceIndex = globalTransforms.index(wfiat);
538 
539  const transformer& transform = globalTransforms.transform
540  (
541  globalTransforms.transformIndex(wfiat)
542  );
543 
544  const face& f = mesh_.faces()[wallFaceIndex];
545 
546  label patchi = mesh_.boundaryMesh().patchID()
547  [
548  wallFaceIndex - mesh_.nInternalFaces()
549  ];
550 
551  referredWallFaces_[rWFI] = referredWallFace
552  (
553  face(identity(f.size())),
554  transform.invTransformPosition(f.points(mesh_.points())),
555  patchi
556  );
557  }
558 
559  wallFaceMap().distribute(referredWallFaces_);
560 
561  if (writeCloud_)
562  {
563  writeReferredWallFaces();
564  }
565 
566  // Direct interaction list and direct wall faces
567 
568  Info<< " Building direct interaction lists" << endl;
569 
570  indexedOctree<treeDataFace> wallFacesTree
571  (
572  treeDataFace(true, mesh_, localWallFaces),
573  procBbRndExt,
574  8, // maxLevel,
575  10, // leafSize,
576  100.0
577  );
578 
579  dil_.setSize(mesh_.nCells());
580 
581  dwfil_.setSize(mesh_.nCells());
582 
583  forAll(cellBbs, celli)
584  {
585  const treeBoundBox& cellBb = cellBbs[celli];
586 
587  treeBoundBox extendedBb
588  (
589  cellBb.min() - interactionVec,
590  cellBb.max() + interactionVec
591  );
592 
593  // Find all cells intersecting extendedBb
594  labelList interactingElems(allCellsTree.findBox(extendedBb));
595 
596  // Reserve space to avoid multiple resizing
597  DynamicList<label> cellDIL(interactingElems.size());
598 
599  forAll(interactingElems, i)
600  {
601  label elemI = interactingElems[i];
602 
603  label c = allCellsTree.shapes().cellLabels()[elemI];
604 
605  // Here, a more detailed geometric test could be applied,
606  // i.e. a more accurate bounding volume like a OBB or
607  // convex hull, or an exact geometrical test.
608 
609  // The higher index cell is added to the lower index
610  // cell's DIL. A cell is not added to its own DIL.
611  if (c > celli)
612  {
613  cellDIL.append(c);
614  }
615  }
616 
617  dil_[celli].transfer(cellDIL);
618 
619  // Find all wall faces intersecting extendedBb
620  interactingElems = wallFacesTree.findBox(extendedBb);
621 
622  dwfil_[celli].setSize(interactingElems.size(), -1);
623 
624  forAll(interactingElems, i)
625  {
626  label elemI = interactingElems[i];
627 
628  label f = wallFacesTree.shapes().faceLabels()[elemI];
629 
630  dwfil_[celli][i] = f;
631  }
632  }
633 }
634 
635 
636 template<class ParticleType>
638 (
639  const treeBoundBox& procBb,
640  const List<treeBoundBox>& allExtendedProcBbs,
641  const globalIndexAndTransform& globalTransforms,
642  List<treeBoundBox>& extendedProcBbsInRange,
643  List<label>& extendedProcBbsTransformIndex,
644  List<label>& extendedProcBbsOrigProc
645 )
646 {
647  extendedProcBbsInRange.setSize(0);
648  extendedProcBbsTransformIndex.setSize(0);
649  extendedProcBbsOrigProc.setSize(0);
650 
651  DynamicList<treeBoundBox> tmpExtendedProcBbsInRange;
652  DynamicList<label> tmpExtendedProcBbsTransformIndex;
653  DynamicList<label> tmpExtendedProcBbsOrigProc;
654 
655  label nTrans = globalTransforms.nIndependentTransforms();
656 
657  forAll(allExtendedProcBbs, proci)
658  {
659  List<label> permutationIndices(nTrans, 0);
660 
661  if (nTrans == 0 && proci != Pstream::myProcNo())
662  {
663  treeBoundBox extendedReferredProcBb = allExtendedProcBbs[proci];
664 
665  if (procBb.overlaps(extendedReferredProcBb))
666  {
667  tmpExtendedProcBbsInRange.append(extendedReferredProcBb);
668 
669  // Dummy index, there are no transforms, so there will
670  // be no resultant transform when this is decoded.
671  tmpExtendedProcBbsTransformIndex.append(0);
672 
673  tmpExtendedProcBbsOrigProc.append(proci);
674  }
675  }
676  else if (nTrans == 3)
677  {
678  label& i = permutationIndices[0];
679  label& j = permutationIndices[1];
680  label& k = permutationIndices[2];
681 
682  for (i = -1; i <= 1; i++)
683  {
684  for (j = -1; j <= 1; j++)
685  {
686  for (k = -1; k <= 1; k++)
687  {
688  if
689  (
690  i == 0
691  && j == 0
692  && k == 0
693  && proci == Pstream::myProcNo()
694  )
695  {
696  // Skip this processor's extended boundBox
697  // when it has no transformation
698  continue;
699  }
700 
701  label transI = globalTransforms.encodeTransformIndex
702  (
703  permutationIndices
704  );
705 
706  const transformer& transform =
707  globalTransforms.transform(transI);
708 
709  treeBoundBox extendedReferredProcBb
710  (
711  transform.transformPosition
712  (
713  allExtendedProcBbs[proci].points()
714  )()
715  );
716 
717  if (procBb.overlaps(extendedReferredProcBb))
718  {
719  tmpExtendedProcBbsInRange.append
720  (
721  extendedReferredProcBb
722  );
723 
724  tmpExtendedProcBbsTransformIndex.append(transI);
725 
726  tmpExtendedProcBbsOrigProc.append(proci);
727  }
728  }
729  }
730  }
731  }
732  else if (nTrans == 2)
733  {
734  label& i = permutationIndices[0];
735  label& j = permutationIndices[1];
736 
737  for (i = -1; i <= 1; i++)
738  {
739  for (j = -1; j <= 1; j++)
740  {
741  if (i == 0 && j == 0 && proci == Pstream::myProcNo())
742  {
743  // Skip this processor's extended boundBox
744  // when it has no transformation
745  continue;
746  }
747 
748  label transI = globalTransforms.encodeTransformIndex
749  (
750  permutationIndices
751  );
752 
753  const transformer& transform =
754  globalTransforms.transform(transI);
755 
756  treeBoundBox extendedReferredProcBb
757  (
758  transform.transformPosition
759  (
760  allExtendedProcBbs[proci].points()
761  )()
762  );
763 
764  if (procBb.overlaps(extendedReferredProcBb))
765  {
766  tmpExtendedProcBbsInRange.append
767  (
768  extendedReferredProcBb
769  );
770 
771  tmpExtendedProcBbsTransformIndex.append(transI);
772 
773  tmpExtendedProcBbsOrigProc.append(proci);
774  }
775  }
776  }
777  }
778  else if (nTrans == 1)
779  {
780  label& i = permutationIndices[0];
781 
782  for (i = -1; i <= 1; i++)
783  {
784  if (i == 0 && proci == Pstream::myProcNo())
785  {
786  // Skip this processor's extended boundBox when it
787  // has no transformation
788  continue;
789  }
790 
791  label transI = globalTransforms.encodeTransformIndex
792  (
793  permutationIndices
794  );
795 
796  const transformer& transform =
797  globalTransforms.transform(transI);
798 
799  treeBoundBox extendedReferredProcBb
800  (
801  transform.transformPosition
802  (
803  allExtendedProcBbs[proci].points()
804  )()
805  );
806 
807  if (procBb.overlaps(extendedReferredProcBb))
808  {
809  tmpExtendedProcBbsInRange.append
810  (
811  extendedReferredProcBb
812  );
813 
814  tmpExtendedProcBbsTransformIndex.append(transI);
815 
816  tmpExtendedProcBbsOrigProc.append(proci);
817  }
818  }
819  }
820  }
821 
822  extendedProcBbsInRange = move(tmpExtendedProcBbsInRange);
823  extendedProcBbsTransformIndex = move(tmpExtendedProcBbsTransformIndex);
824  extendedProcBbsOrigProc = move(tmpExtendedProcBbsOrigProc);
825 }
826 
827 
828 template<class ParticleType>
830 (
831  autoPtr<distributionMap>& mapPtr,
832  const List<label>& toProc
833 )
834 {
835  // Determine send map
836  // ~~~~~~~~~~~~~~~~~~
837 
838  // 1. Count
839  labelList nSend(Pstream::nProcs(), 0);
840 
841  forAll(toProc, i)
842  {
843  label proci = toProc[i];
844 
845  nSend[proci]++;
846  }
847 
848  // 2. Size sendMap
849  labelListList sendMap(Pstream::nProcs());
850 
851  forAll(nSend, proci)
852  {
853  sendMap[proci].setSize(nSend[proci]);
854 
855  nSend[proci] = 0;
856  }
857 
858  // 3. Fill sendMap
859  forAll(toProc, i)
860  {
861  label proci = toProc[i];
862 
863  sendMap[proci][nSend[proci]++] = i;
864  }
865 
866  // 4. Send over how many I need to receive
867  labelList recvSizes;
868  Pstream::exchangeSizes(sendMap, recvSizes);
869 
870 
871  // Determine receive map
872  // ~~~~~~~~~~~~~~~~~~~~~
873 
874  labelListList constructMap(Pstream::nProcs());
875 
876  // Local transfers first
877  constructMap[Pstream::myProcNo()] = identity
878  (
879  sendMap[Pstream::myProcNo()].size()
880  );
881 
882  label constructSize = constructMap[Pstream::myProcNo()].size();
883 
884  forAll(constructMap, proci)
885  {
886  if (proci != Pstream::myProcNo())
887  {
888  label nRecv = recvSizes[proci];
889 
890  constructMap[proci].setSize(nRecv);
891 
892  for (label i = 0; i < nRecv; i++)
893  {
894  constructMap[proci][i] = constructSize++;
895  }
896  }
897  }
898 
899  mapPtr.reset
900  (
901  new distributionMap
902  (
903  constructSize,
904  move(sendMap),
905  move(constructMap)
906  )
907  );
908 }
909 
910 
911 template<class ParticleType>
913 (
914  const List<DynamicList<ParticleType*>>& cellOccupancy
915 )
916 {
917  const globalIndexAndTransform& globalTransforms =
918  mesh_.globalData().globalTransforms();
919 
920  referredParticles_.setSize(cellIndexAndTransformToDistribute_.size());
921 
922  // Clear all existing referred particles
923 
924  forAll(referredParticles_, i)
925  {
926  referredParticles_[i].clear();
927  }
928 
929  // Clear all particles that may have been populated into the cloud
930  cloud_.clear();
931 
932  forAll(cellIndexAndTransformToDistribute_, i)
933  {
934  const labelPair ciat = cellIndexAndTransformToDistribute_[i];
935 
936  label cellIndex = globalTransforms.index(ciat);
937 
938  List<ParticleType*> realParticles = cellOccupancy[cellIndex];
939 
940  IDLList<ParticleType>& particlesToRefer = referredParticles_[i];
941 
942  forAll(realParticles, rM)
943  {
944  const ParticleType& particle = *realParticles[rM];
945 
946  particlesToRefer.append(particle.clone().ptr());
947 
948  prepareParticleToBeReferred(particlesToRefer.last(), ciat);
949  }
950  }
951 }
952 
953 
954 template<class ParticleType>
956 (
957  ParticleType* particle,
958  labelPair ciat
959 )
960 {
961  const globalIndexAndTransform& globalTransforms =
962  mesh_.globalData().globalTransforms();
963 
964  const transformer& transform = globalTransforms.transform
965  (
966  globalTransforms.transformIndex(ciat)
967  );
968 
969  particle->prepareForInteractionListReferral(transform);
970 }
971 
972 
973 template<class ParticleType>
975 {
976  if (writeCloud_)
977  {
978  forAll(referredParticles_, refCelli)
979  {
980  const IDLList<ParticleType>& refCell =
981  referredParticles_[refCelli];
982 
983  forAllConstIter(typename IDLList<ParticleType>, refCell, iter)
984  {
985  cloud_.addParticle
986  (
987  static_cast<ParticleType*>(iter().clone().ptr())
988  );
989  }
990  }
991  }
992 }
993 
994 
995 template<class ParticleType>
997 {
998  const globalIndexAndTransform& globalTransforms =
999  mesh_.globalData().globalTransforms();
1000 
1001  referredWallData_.setSize
1002  (
1003  wallFaceIndexAndTransformToDistribute_.size()
1004  );
1005 
1006  const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
1007 
1008  forAll(referredWallData_, rWVI)
1009  {
1010  const labelPair& wfiat = wallFaceIndexAndTransformToDistribute_[rWVI];
1011 
1012  label wallFaceIndex = globalTransforms.index(wfiat);
1013 
1014  const transformer& transform = globalTransforms.transform
1015  (
1016  globalTransforms.transformIndex(wfiat)
1017  );
1018 
1019  label patchi = mesh_.boundaryMesh().patchID()
1020  [
1021  wallFaceIndex - mesh_.nInternalFaces()
1022  ];
1023 
1024  label patchFacei =
1025  wallFaceIndex
1026  - mesh_.boundaryMesh()[patchi].start();
1027 
1028  // Need to transform velocity when tensor transforms are
1029  // supported
1030  referredWallData_[rWVI] = U.boundaryField()[patchi][patchFacei];
1031 
1032  if (transform.transforms())
1033  {
1034  referredWallData_[rWVI] =
1035  transform.invTransform(referredWallData_[rWVI]);
1036  }
1037  }
1038 }
1039 
1040 
1041 template<class ParticleType>
1043 {
1044  if (referredWallFaces_.empty())
1045  {
1046  return;
1047  }
1048 
1049  fileName objDir = mesh_.time().timePath()/cloud::prefix;
1050 
1051  mkDir(objDir);
1052 
1053  fileName objFileName = "referredWallFaces.obj";
1054 
1055  OFstream str(objDir/objFileName);
1056 
1057  Info<< " Writing "
1058  << mesh_.time().timeName()/cloud::prefix/objFileName
1059  << endl;
1060 
1061  label offset = 1;
1062 
1063  forAll(referredWallFaces_, rWFI)
1064  {
1065  const referredWallFace& rwf = referredWallFaces_[rWFI];
1066 
1067  forAll(rwf, fPtI)
1068  {
1069  meshTools::writeOBJ(str, rwf.points()[rwf[fPtI]]);
1070  }
1071 
1072  str<< 'f';
1073 
1074  forAll(rwf, fPtI)
1075  {
1076  str<< ' ' << fPtI + offset;
1077  }
1078 
1079  str<< nl;
1080 
1081  offset += rwf.size();
1082  }
1083 }
1084 
1085 
1086 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1087 
1088 template<class ParticleType>
1090 :
1091  mesh_(mesh),
1092  cloud_(mesh_, "nullptr_Cloud", IDLList<ParticleType>()),
1093  writeCloud_(false),
1094  cellMapPtr_(),
1095  wallFaceMapPtr_(),
1096  maxDistance_(0.0),
1097  dil_(),
1098  dwfil_(),
1099  ril_(),
1100  rilInverse_(),
1101  cellIndexAndTransformToDistribute_(),
1102  wallFaceIndexAndTransformToDistribute_(),
1103  referredWallFaces_(),
1104  UName_("unknown_U"),
1105  referredWallData_(),
1106  referredParticles_()
1107 {}
1108 
1109 
1110 template<class ParticleType>
1113  const polyMesh& mesh,
1114  scalar maxDistance,
1115  Switch writeCloud,
1116  const word& UName
1117 )
1118 :
1119  mesh_(mesh),
1120  cloud_(mesh_, "referredParticleCloud", IDLList<ParticleType>()),
1121  writeCloud_(writeCloud),
1122  cellMapPtr_(),
1123  wallFaceMapPtr_(),
1124  maxDistance_(maxDistance),
1125  dil_(),
1126  dwfil_(),
1127  ril_(),
1128  rilInverse_(),
1129  cellIndexAndTransformToDistribute_(),
1130  wallFaceIndexAndTransformToDistribute_(),
1131  referredWallFaces_(),
1132  UName_(UName),
1133  referredWallData_(),
1134  referredParticles_()
1135 {
1136  buildInteractionLists();
1137 }
1138 
1139 
1140 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1141 
1142 template<class ParticleType>
1144 {}
1145 
1146 
1147 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
1148 
1149 template<class ParticleType>
1153  PstreamBuffers& pBufs
1154 )
1155 {
1156  if (mesh_.changing())
1157  {
1159  << "Mesh changing, rebuilding InteractionLists form scratch."
1160  << endl;
1161 
1162  buildInteractionLists();
1163  }
1164 
1165  prepareWallDataToRefer();
1166 
1167  prepareParticlesToRefer(cellOccupancy);
1168 
1169  for (label domain = 0; domain < Pstream::nProcs(); domain++)
1170  {
1171  const labelList& subMap = cellMap().subMap()[domain];
1172 
1173  if (subMap.size())
1174  {
1175  UOPstream toDomain(domain, pBufs);
1176 
1177  UIndirectList<IDLList<ParticleType>> subMappedParticles
1178  (
1179  referredParticles_,
1180  subMap
1181  );
1182 
1183  forAll(subMappedParticles, i)
1184  {
1185  toDomain << subMappedParticles[i];
1186  }
1187  }
1188  }
1189 
1190  // Using the distributionMap to start sending and receiving the
1191  // buffer but not block, i.e. it is calling
1192  // pBufs.finishedSends(false);
1193  wallFaceMap().send(pBufs, referredWallData_);
1194 }
1195 
1196 
1197 template<class ParticleType>
1200  PstreamBuffers& pBufs,
1201  const label startOfRequests
1202 )
1203 {
1204  Pstream::waitRequests(startOfRequests);
1205 
1206  referredParticles_.setSize(cellMap().constructSize());
1207 
1208  for (label domain = 0; domain < Pstream::nProcs(); domain++)
1209  {
1210  const labelList& constructMap = cellMap().constructMap()[domain];
1211 
1212  if (constructMap.size())
1213  {
1214  UIPstream str(domain, pBufs);
1215 
1216  forAll(constructMap, i)
1217  {
1218  referredParticles_[constructMap[i]] = IDLList<ParticleType>
1219  (
1220  str,
1221  typename ParticleType::iNew(mesh_)
1222  );
1223  }
1224  }
1225  }
1226 
1227  forAll(referredParticles_, refCelli)
1228  {
1229  IDLList<ParticleType>& refCell = referredParticles_[refCelli];
1230  forAllIter(typename IDLList<ParticleType>, refCell, iter)
1231  {
1232  iter().correctAfterInteractionListReferral(ril_[refCelli][0]);
1233  }
1234  }
1235 
1236  fillReferredParticleCloud();
1237 
1238  wallFaceMap().receive(pBufs, referredWallData_);
1239 }
1240 
1241 
1242 // ************************************************************************* //
Template class for intrusive linked lists.
Definition: ILList.H:50
void inplaceSubset(const UList< T > &select, const T &value, ListType &)
Inplace extract elements of List when select is a certain value.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
List< treeBoundBox > treeBoundBoxList
List of bounding boxes.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Builds direct interaction list, specifying which local (real) cells are potentially in range of each ...
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none/any.
Definition: Switch.H:60
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:59
const dimensionedScalar c
Speed of light in a vacuum.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
Input inter-processor communications stream operating on external buffer.
Definition: UIPstream.H:53
T clone(const T &t)
Definition: List.H:54
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
A class for handling words, derived from string.
Definition: word.H:59
mkDir(pdfPath)
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
List< label > labelList
A List of labels.
Definition: labelList.H:56
Output inter-processor communications stream operating on external buffer.
Definition: UOPstream.H:54
void receiveReferredData(PstreamBuffers &pBufs, const label startReq=0)
Receive referred data.
static const char nl
Definition: Ostream.H:260
labelList f(nPoints)
Buffers for inter-processor communications streams (UOPstream, UIPstream).
void setSize(const label)
Reset size of List.
Definition: List.C:281
label patchi
const List< DynamicList< molecule * > > & cellOccupancy
#define WarningInFunction
Report a warning using Foam::Warning.
A List with indirect addressing.
Definition: fvMatrix.H:106
void transform(const fvPatch &patch, const label patchFacei, const transformer &transform, TrackingData &td)
Transform across an interface.
messageStream Info
InteractionLists(const polyMesh &mesh)
Construct null from mesh.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105
void sendReferredData(const List< DynamicList< ParticleType *>> &cellOccupancy, PstreamBuffers &pBufs)
Prepare and send referred particles and wall data,.