meshToMeshParallelOps.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) 2012-2017 OpenFOAM Foundation
9  Copyright (C) 2015-2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "meshToMesh.H"
30 #include "OFstream.H"
31 #include "Time.H"
32 #include "globalIndex.H"
33 #include "mergePoints.H"
34 #include "processorPolyPatch.H"
35 #include "SubField.H"
36 #include "AABBTree.H"
37 #include "cellBox.H"
38 
39 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40 
41 Foam::label Foam::meshToMesh::calcDistribution
42 (
43  const polyMesh& src,
44  const polyMesh& tgt
45 ) const
46 {
47  label proci = 0;
48 
49  if (Pstream::parRun())
50  {
51  List<label> cellsPresentOnProc(Pstream::nProcs(), Zero);
52  if ((src.nCells() > 0) || (tgt.nCells() > 0))
53  {
54  cellsPresentOnProc[Pstream::myProcNo()] = 1;
55  }
56  else
57  {
58  cellsPresentOnProc[Pstream::myProcNo()] = 0;
59  }
60 
61  Pstream::gatherList(cellsPresentOnProc);
62  Pstream::scatterList(cellsPresentOnProc);
63 
64  label nHaveCells = sum(cellsPresentOnProc);
65 
66  if (nHaveCells > 1)
67  {
68  proci = -1;
69 
70  DebugInFunction
71  << "Meshes split across multiple processors" << endl;
72  }
73  else if (nHaveCells == 1)
74  {
75  proci = cellsPresentOnProc.find(1);
76 
77  DebugInFunction
78  << "Meshes local to processor" << proci << endl;
79  }
80  }
81 
82  return proci;
83 }
84 
85 
86 Foam::label Foam::meshToMesh::calcOverlappingProcs
87 (
88  const List<treeBoundBoxList>& procBb,
89  const boundBox& bb,
90  boolList& overlaps
91 ) const
92 {
93  overlaps = false;
94 
95  label nOverlaps = 0;
96 
97  forAll(procBb, proci)
98  {
99  const treeBoundBoxList& bbp = procBb[proci];
100 
101  for (const treeBoundBox& b : bbp)
102  {
103  if (b.overlaps(bb))
104  {
105  overlaps[proci] = true;
106  ++nOverlaps;
107  break;
108  }
109  }
110  }
111 
112  return nOverlaps;
113 }
114 
115 
116 Foam::autoPtr<Foam::mapDistribute> Foam::meshToMesh::calcProcMap
117 (
118  const polyMesh& src,
119  const polyMesh& tgt
120 ) const
121 {
122  switch (procMapMethod_)
123  {
124  case procMapMethod::pmLOD:
125  {
126  Info<< "meshToMesh: Using processorLOD method" << endl;
127 
128  // Create processor map of overlapping faces. This map gets
129  // (possibly remote) cells from the tgt mesh such that they
130  // (together) cover all of the src mesh
131  const label nGlobalSrcCells = src.globalData().nTotalCells();
132  const label cellsPerBox = max(1, 0.001*nGlobalSrcCells);
133  typename processorLODs::cellBox boxLOD
134  (
135  src.cells(),
136  src.faces(),
137  src.points(),
138  tgt.cells(),
139  tgt.faces(),
140  tgt.points(),
141  cellsPerBox,
142  src.nCells()
143  );
144 
145  return boxLOD.map();
146  break;
147  }
148  default:
149  {
150  Info<< "meshToMesh: Using AABBTree method" << endl;
151 
152  // get decomposition of cells on src mesh
153  List<treeBoundBoxList> procBb(Pstream::nProcs());
154 
155  if (src.nCells() > 0)
156  {
157  procBb[Pstream::myProcNo()] = AABBTree<labelList>
158  (
159  src.cellPoints(),
160  src.points(),
161  false
162  ).boundBoxes();
163  }
164  else
165  {
166  procBb[Pstream::myProcNo()] = treeBoundBoxList();
167  }
168 
169 
170  Pstream::gatherList(procBb);
171  Pstream::scatterList(procBb);
172 
173 
174  if (debug)
175  {
176  InfoInFunction
177  << "Determining extent of src mesh per processor:" << nl
178  << "\tproc\tbb" << endl;
179  forAll(procBb, proci)
180  {
181  Info<< '\t' << proci << '\t' << procBb[proci] << endl;
182  }
183  }
184 
185 
186  // determine which cells of tgt mesh overlaps src mesh per proc
187  const cellList& cells = tgt.cells();
188  const faceList& faces = tgt.faces();
189  const pointField& points = tgt.points();
190 
191  labelListList sendMap;
192 
193  {
194  // per processor indices into all segments to send
195  List<DynamicList<label>> dynSendMap(Pstream::nProcs());
196  label iniSize = floor(tgt.nCells()/Pstream::nProcs());
197 
198  forAll(dynSendMap, proci)
199  {
200  dynSendMap[proci].setCapacity(iniSize);
201  }
202 
203  // work array - whether src processor bb overlaps the tgt cell
204  // bounds
205  boolList procBbOverlaps(Pstream::nProcs());
206  forAll(cells, celli)
207  {
208  const cell& c = cells[celli];
209 
210  // determine bounding box of tgt cell
211  boundBox cellBb(boundBox::invertedBox);
212  forAll(c, facei)
213  {
214  const face& f = faces[c[facei]];
215  forAll(f, fp)
216  {
217  cellBb.add(points, f);
218  }
219  }
220 
221  // find the overlapping tgt cells on each src processor
222  (void)calcOverlappingProcs(procBb, cellBb, procBbOverlaps);
223 
224  forAll(procBbOverlaps, proci)
225  {
226  if (procBbOverlaps[proci])
227  {
228  dynSendMap[proci].append(celli);
229  }
230  }
231  }
232 
233  // convert dynamicList to labelList
234  sendMap.setSize(Pstream::nProcs());
235  forAll(sendMap, proci)
236  {
237  sendMap[proci].transfer(dynSendMap[proci]);
238  }
239  }
240 
241  // debug printing
242  if (debug)
243  {
244  Pout<< "Of my " << cells.size()
245  << " target cells I need to send to:" << nl
246  << "\tproc\tcells" << endl;
247  forAll(sendMap, proci)
248  {
249  Pout<< '\t' << proci << '\t' << sendMap[proci].size()
250  << endl;
251  }
252  }
253 
254 
255  // send over how many tgt cells I need to receive from each
256  // processor
257  labelListList sendSizes(Pstream::nProcs());
258  sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
259  forAll(sendMap, proci)
260  {
261  sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
262  }
263  Pstream::gatherList(sendSizes);
264  Pstream::scatterList(sendSizes);
265 
266 
267  // determine order of receiving
268  labelListList constructMap(Pstream::nProcs());
269 
270  label segmentI = 0;
271  forAll(constructMap, proci)
272  {
273  // what I need to receive is what other processor is sending
274  // to me
275  label nRecv = sendSizes[proci][Pstream::myProcNo()];
276  constructMap[proci].setSize(nRecv);
277 
278  for (label i = 0; i < nRecv; i++)
279  {
280  constructMap[proci][i] = segmentI++;
281  }
282  }
283 
284  return autoPtr<mapDistribute>::New
285  (
286  segmentI, // size after construction
287  std::move(sendMap),
288  std::move(constructMap)
289  );
290 
291  break;
292  }
293  }
294 }
295 
296 
297 void Foam::meshToMesh::distributeCells
298 (
299  const mapDistribute& map,
300  const polyMesh& tgtMesh,
301  const globalIndex& globalI,
302  List<pointField>& points,
303  List<label>& nInternalFaces,
304  List<faceList>& faces,
305  List<labelList>& faceOwner,
306  List<labelList>& faceNeighbour,
307  List<labelList>& cellIDs,
308  List<labelList>& nbrProcIDs,
309  List<labelList>& procLocalFaceIDs
310 ) const
311 {
312  PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
313 
314  points.setSize(Pstream::nProcs());
315  nInternalFaces.setSize(Pstream::nProcs(), 0);
316  faces.setSize(Pstream::nProcs());
317  faceOwner.setSize(Pstream::nProcs());
318  faceNeighbour.setSize(Pstream::nProcs());
319  cellIDs.setSize(Pstream::nProcs());
320 
321  nbrProcIDs.setSize(Pstream::nProcs());;
322  procLocalFaceIDs.setSize(Pstream::nProcs());;
323 
324 
325  for (const int domain : Pstream::allProcs())
326  {
327  const labelList& sendElems = map.subMap()[domain];
328 
329  if (sendElems.size())
330  {
331  // reverse cell map
332  labelList reverseCellMap(tgtMesh.nCells(), -1);
333  forAll(sendElems, subCelli)
334  {
335  reverseCellMap[sendElems[subCelli]] = subCelli;
336  }
337 
338  DynamicList<face> subFaces(tgtMesh.nFaces());
339  DynamicList<label> subFaceOwner(tgtMesh.nFaces());
340  DynamicList<label> subFaceNeighbour(tgtMesh.nFaces());
341 
342  DynamicList<label> subNbrProcIDs(tgtMesh.nFaces());
343  DynamicList<label> subProcLocalFaceIDs(tgtMesh.nFaces());
344 
345  label nInternal = 0;
346 
347  // internal faces
348  forAll(tgtMesh.faceNeighbour(), facei)
349  {
350  label own = tgtMesh.faceOwner()[facei];
351  label nbr = tgtMesh.faceNeighbour()[facei];
352  label subOwn = reverseCellMap[own];
353  label subNbr = reverseCellMap[nbr];
354 
355  if (subOwn != -1 && subNbr != -1)
356  {
357  nInternal++;
358 
359  if (subOwn < subNbr)
360  {
361  subFaces.append(tgtMesh.faces()[facei]);
362  subFaceOwner.append(subOwn);
363  subFaceNeighbour.append(subNbr);
364  subNbrProcIDs.append(-1);
365  subProcLocalFaceIDs.append(-1);
366  }
367  else
368  {
369  subFaces.append(tgtMesh.faces()[facei].reverseFace());
370  subFaceOwner.append(subNbr);
371  subFaceNeighbour.append(subOwn);
372  subNbrProcIDs.append(-1);
373  subProcLocalFaceIDs.append(-1);
374  }
375  }
376  }
377 
378  // boundary faces for new region
379  forAll(tgtMesh.faceNeighbour(), facei)
380  {
381  label own = tgtMesh.faceOwner()[facei];
382  label nbr = tgtMesh.faceNeighbour()[facei];
383  label subOwn = reverseCellMap[own];
384  label subNbr = reverseCellMap[nbr];
385 
386  if (subOwn != -1 && subNbr == -1)
387  {
388  subFaces.append(tgtMesh.faces()[facei]);
389  subFaceOwner.append(subOwn);
390  subFaceNeighbour.append(subNbr);
391  subNbrProcIDs.append(-1);
392  subProcLocalFaceIDs.append(-1);
393  }
394  else if (subOwn == -1 && subNbr != -1)
395  {
396  subFaces.append(tgtMesh.faces()[facei].reverseFace());
397  subFaceOwner.append(subNbr);
398  subFaceNeighbour.append(subOwn);
399  subNbrProcIDs.append(-1);
400  subProcLocalFaceIDs.append(-1);
401  }
402  }
403 
404  // boundary faces of existing region
405  forAll(tgtMesh.boundaryMesh(), patchi)
406  {
407  const polyPatch& pp = tgtMesh.boundaryMesh()[patchi];
408 
409  label nbrProci = -1;
410 
411  // store info for faces on processor patches
412  if (isA<processorPolyPatch>(pp))
413  {
414  const processorPolyPatch& ppp =
415  dynamic_cast<const processorPolyPatch&>(pp);
416 
417  nbrProci = ppp.neighbProcNo();
418  }
419 
420  forAll(pp, i)
421  {
422  label facei = pp.start() + i;
423  label own = tgtMesh.faceOwner()[facei];
424 
425  if (reverseCellMap[own] != -1)
426  {
427  subFaces.append(tgtMesh.faces()[facei]);
428  subFaceOwner.append(reverseCellMap[own]);
429  subFaceNeighbour.append(-1);
430  subNbrProcIDs.append(nbrProci);
431  subProcLocalFaceIDs.append(i);
432  }
433  }
434  }
435 
436  // reverse point map
437  labelList reversePointMap(tgtMesh.nPoints(), -1);
438  DynamicList<point> subPoints(tgtMesh.nPoints());
439  forAll(subFaces, subFacei)
440  {
441  face& f = subFaces[subFacei];
442  forAll(f, fp)
443  {
444  label pointi = f[fp];
445  if (reversePointMap[pointi] == -1)
446  {
447  reversePointMap[pointi] = subPoints.size();
448  subPoints.append(tgtMesh.points()[pointi]);
449  }
450 
451  f[fp] = reversePointMap[pointi];
452  }
453  }
454 
455  // tgt cells into global numbering
456  labelList globalElems(globalI.toGlobal(sendElems));
457 
458  if (debug > 1)
459  {
460  forAll(sendElems, i)
461  {
462  Pout<< "tgtProc:" << Pstream::myProcNo()
463  << " sending tgt cell " << sendElems[i]
464  << "[" << globalElems[i] << "]"
465  << " to srcProc " << domain << endl;
466  }
467  }
468 
469  // pass data
470  if (domain == Pstream::myProcNo())
471  {
472  // allocate my own data
473  points[Pstream::myProcNo()] = subPoints;
474  nInternalFaces[Pstream::myProcNo()] = nInternal;
475  faces[Pstream::myProcNo()] = subFaces;
476  faceOwner[Pstream::myProcNo()] = subFaceOwner;
477  faceNeighbour[Pstream::myProcNo()] = subFaceNeighbour;
478  cellIDs[Pstream::myProcNo()] = globalElems;
479  nbrProcIDs[Pstream::myProcNo()] = subNbrProcIDs;
480  procLocalFaceIDs[Pstream::myProcNo()] = subProcLocalFaceIDs;
481  }
482  else
483  {
484  // send data to other processor domains
485  UOPstream toDomain(domain, pBufs);
486 
487  toDomain
488  << subPoints
489  << nInternal
490  << subFaces
491  << subFaceOwner
492  << subFaceNeighbour
493  << globalElems
494  << subNbrProcIDs
495  << subProcLocalFaceIDs;
496  }
497  }
498  }
499 
500  // Start receiving
501  pBufs.finishedSends();
502 
503  // Consume
504  for (const int domain : Pstream::allProcs())
505  {
506  const labelList& recvElems = map.constructMap()[domain];
507 
508  if (domain != Pstream::myProcNo() && recvElems.size())
509  {
510  UIPstream str(domain, pBufs);
511 
512  str >> points[domain]
513  >> nInternalFaces[domain]
514  >> faces[domain]
515  >> faceOwner[domain]
516  >> faceNeighbour[domain]
517  >> cellIDs[domain]
518  >> nbrProcIDs[domain]
519  >> procLocalFaceIDs[domain];
520  }
521 
522  if (debug)
523  {
524  Pout<< "Target mesh send sizes[" << domain << "]"
525  << ": points="<< points[domain].size()
526  << ", faces=" << faces[domain].size()
527  << ", nInternalFaces=" << nInternalFaces[domain]
528  << ", faceOwn=" << faceOwner[domain].size()
529  << ", faceNbr=" << faceNeighbour[domain].size()
530  << ", cellIDs=" << cellIDs[domain].size() << endl;
531  }
532  }
533 }
534 
535 
536 void Foam::meshToMesh::distributeAndMergeCells
537 (
538  const mapDistribute& map,
539  const polyMesh& tgt,
540  const globalIndex& globalI,
541  pointField& tgtPoints,
542  faceList& tgtFaces,
543  labelList& tgtFaceOwners,
544  labelList& tgtFaceNeighbours,
545  labelList& tgtCellIDs
546 ) const
547 {
548  // Exchange per-processor data
549  List<pointField> allPoints;
550  List<label> allNInternalFaces;
551  List<faceList> allFaces;
552  List<labelList> allFaceOwners;
553  List<labelList> allFaceNeighbours;
554  List<labelList> allTgtCellIDs;
555 
556  // Per target mesh face the neighbouring proc and index in
557  // processor patch (all -1 for normal boundary face)
558  List<labelList> allNbrProcIDs;
559  List<labelList> allProcLocalFaceIDs;
560 
561  distributeCells
562  (
563  map,
564  tgt,
565  globalI,
566  allPoints,
567  allNInternalFaces,
568  allFaces,
569  allFaceOwners,
570  allFaceNeighbours,
571  allTgtCellIDs,
572  allNbrProcIDs,
573  allProcLocalFaceIDs
574  );
575 
576  // Convert lists into format that can be used to generate a valid polyMesh
577  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
578  //
579  // Points and cells are collected into single flat lists:
580  // - i.e. proc0, proc1 ... procN
581  //
582  // Faces need to be sorted after collection to that internal faces are
583  // contiguous, followed by all boundary faces
584  //
585  // Processor patch faces between included cells on neighbouring processors
586  // are converted into internal faces
587  //
588  // Face list structure:
589  // - Per processor:
590  // - internal faces
591  // - processor faces that have been converted into internal faces
592  // - Followed by all boundary faces
593  // - from 'normal' boundary faces
594  // - from singularly-sided processor patch faces
595 
596 
597  // Number of internal+coupled faces
598  labelList allNIntCoupledFaces(allNInternalFaces);
599 
600  // Starting offset for points
601  label nPoints = 0;
602  labelList pointOffset(Pstream::nProcs(), Zero);
603  forAll(allPoints, proci)
604  {
605  pointOffset[proci] = nPoints;
606  nPoints += allPoints[proci].size();
607  }
608 
609  // Starting offset for cells
610  label nCells = 0;
611  labelList cellOffset(Pstream::nProcs(), Zero);
612  forAll(allTgtCellIDs, proci)
613  {
614  cellOffset[proci] = nCells;
615  nCells += allTgtCellIDs[proci].size();
616  }
617 
618  // Count any coupled faces
619  typedef FixedList<label, 3> label3;
620  typedef HashTable<label, label3, label3::Hash<>> procCoupleInfo;
621  procCoupleInfo procFaceToGlobalCell;
622 
623  forAll(allNbrProcIDs, proci)
624  {
625  const labelList& nbrProci = allNbrProcIDs[proci];
626  const labelList& localFacei = allProcLocalFaceIDs[proci];
627 
628  forAll(nbrProci, i)
629  {
630  if (nbrProci[i] != -1 && localFacei[i] != -1)
631  {
632  label3 key;
633  key[0] = min(proci, nbrProci[i]);
634  key[1] = max(proci, nbrProci[i]);
635  key[2] = localFacei[i];
636 
637  const auto fnd = procFaceToGlobalCell.cfind(key);
638 
639  if (!fnd.found())
640  {
641  procFaceToGlobalCell.insert(key, -1);
642  }
643  else
644  {
645  if (debug > 1)
646  {
647  Pout<< "Additional internal face between procs:"
648  << key[0] << " and " << key[1]
649  << " across local face " << key[2] << endl;
650  }
651 
652  allNIntCoupledFaces[proci]++;
653  }
654  }
655  }
656  }
657 
658 
659  // Starting offset for internal faces
660  label nIntFaces = 0;
661  label nFacesTotal = 0;
662  labelList internalFaceOffset(Pstream::nProcs(), Zero);
663  forAll(allNIntCoupledFaces, proci)
664  {
665  label nCoupledFaces =
666  allNIntCoupledFaces[proci] - allNInternalFaces[proci];
667 
668  internalFaceOffset[proci] = nIntFaces;
669  nIntFaces += allNIntCoupledFaces[proci];
670  nFacesTotal += allFaceOwners[proci].size() - nCoupledFaces;
671  }
672 
673  tgtPoints.setSize(nPoints);
674  tgtFaces.setSize(nFacesTotal);
675  tgtFaceOwners.setSize(nFacesTotal);
676  tgtFaceNeighbours.setSize(nFacesTotal);
677  tgtCellIDs.setSize(nCells);
678 
679  // Insert points
680  forAll(allPoints, proci)
681  {
682  const pointField& pts = allPoints[proci];
683  SubList<point>(tgtPoints, pts.size(), pointOffset[proci]) = pts;
684  }
685 
686  // Insert cellIDs
687  forAll(allTgtCellIDs, proci)
688  {
689  const labelList& cellIDs = allTgtCellIDs[proci];
690  SubList<label>(tgtCellIDs, cellIDs.size(), cellOffset[proci]) = cellIDs;
691  }
692 
693 
694  // Insert internal faces (from internal faces)
695  forAll(allFaces, proci)
696  {
697  const faceList& fcs = allFaces[proci];
698  const labelList& faceOs = allFaceOwners[proci];
699  const labelList& faceNs = allFaceNeighbours[proci];
700 
701  SubList<face> slice
702  (
703  tgtFaces,
704  allNInternalFaces[proci],
705  internalFaceOffset[proci]
706  );
707  slice = SubList<face>(fcs, allNInternalFaces[proci]);
708  forAll(slice, i)
709  {
710  add(slice[i], pointOffset[proci]);
711  }
712 
713  SubField<label> ownSlice
714  (
715  tgtFaceOwners,
716  allNInternalFaces[proci],
717  internalFaceOffset[proci]
718  );
719  ownSlice = SubField<label>(faceOs, allNInternalFaces[proci]);
720  add(ownSlice, cellOffset[proci]);
721 
722  SubField<label> nbrSlice
723  (
724  tgtFaceNeighbours,
725  allNInternalFaces[proci],
726  internalFaceOffset[proci]
727  );
728  nbrSlice = SubField<label>(faceNs, allNInternalFaces[proci]);
729  add(nbrSlice, cellOffset[proci]);
730 
731  internalFaceOffset[proci] += allNInternalFaces[proci];
732  }
733 
734 
735  // Insert internal faces (from coupled face-pairs)
736  forAll(allNbrProcIDs, proci)
737  {
738  const labelList& nbrProci = allNbrProcIDs[proci];
739  const labelList& localFacei = allProcLocalFaceIDs[proci];
740  const labelList& faceOs = allFaceOwners[proci];
741  const faceList& fcs = allFaces[proci];
742 
743  forAll(nbrProci, i)
744  {
745  if (nbrProci[i] != -1 && localFacei[i] != -1)
746  {
747  label3 key;
748  key[0] = min(proci, nbrProci[i]);
749  key[1] = max(proci, nbrProci[i]);
750  key[2] = localFacei[i];
751 
752  auto fnd = procFaceToGlobalCell.find(key);
753 
754  if (fnd.found())
755  {
756  label tgtFacei = fnd();
757  if (tgtFacei == -1)
758  {
759  // on first visit store the new cell on this side
760  fnd() = cellOffset[proci] + faceOs[i];
761  }
762  else
763  {
764  // get owner and neighbour in new cell numbering
765  label newOwn = cellOffset[proci] + faceOs[i];
766  label newNbr = fnd();
767  label tgtFacei = internalFaceOffset[proci]++;
768 
769  if (debug > 1)
770  {
771  Pout<< " proc " << proci
772  << "\tinserting face:" << tgtFacei
773  << " connection between owner " << newOwn
774  << " and neighbour " << newNbr
775  << endl;
776  }
777 
778  if (newOwn < newNbr)
779  {
780  // we have correct orientation
781  tgtFaces[tgtFacei] = fcs[i];
782  tgtFaceOwners[tgtFacei] = newOwn;
783  tgtFaceNeighbours[tgtFacei] = newNbr;
784  }
785  else
786  {
787  // reverse orientation
788  tgtFaces[tgtFacei] = fcs[i].reverseFace();
789  tgtFaceOwners[tgtFacei] = newNbr;
790  tgtFaceNeighbours[tgtFacei] = newOwn;
791  }
792 
793  add(tgtFaces[tgtFacei], pointOffset[proci]);
794 
795  // mark with unique value
796  fnd() = -2;
797  }
798  }
799  }
800  }
801  }
802 
803 
804  forAll(allNbrProcIDs, proci)
805  {
806  const labelList& nbrProci = allNbrProcIDs[proci];
807  const labelList& localFacei = allProcLocalFaceIDs[proci];
808  const labelList& faceOs = allFaceOwners[proci];
809  const labelList& faceNs = allFaceNeighbours[proci];
810  const faceList& fcs = allFaces[proci];
811 
812  forAll(nbrProci, i)
813  {
814  // coupled boundary face
815  if (nbrProci[i] != -1 && localFacei[i] != -1)
816  {
817  label3 key;
818  key[0] = min(proci, nbrProci[i]);
819  key[1] = max(proci, nbrProci[i]);
820  key[2] = localFacei[i];
821 
822  label tgtFacei = procFaceToGlobalCell[key];
823 
824  if (tgtFacei == -1)
825  {
826  FatalErrorInFunction
827  << "Unvisited " << key
828  << abort(FatalError);
829  }
830  else if (tgtFacei != -2)
831  {
832  label newOwn = cellOffset[proci] + faceOs[i];
833  label tgtFacei = nIntFaces++;
834 
835  if (debug > 1)
836  {
837  Pout<< " proc " << proci
838  << "\tinserting boundary face:" << tgtFacei
839  << " from coupled face " << key
840  << endl;
841  }
842 
843  tgtFaces[tgtFacei] = fcs[i];
844  add(tgtFaces[tgtFacei], pointOffset[proci]);
845 
846  tgtFaceOwners[tgtFacei] = newOwn;
847  tgtFaceNeighbours[tgtFacei] = -1;
848  }
849  }
850  // normal boundary face
851  else
852  {
853  label own = faceOs[i];
854  label nbr = faceNs[i];
855  if ((own != -1) && (nbr == -1))
856  {
857  label newOwn = cellOffset[proci] + faceOs[i];
858  label tgtFacei = nIntFaces++;
859 
860  tgtFaces[tgtFacei] = fcs[i];
861  add(tgtFaces[tgtFacei], pointOffset[proci]);
862 
863  tgtFaceOwners[tgtFacei] = newOwn;
864  tgtFaceNeighbours[tgtFacei] = -1;
865  }
866  }
867  }
868  }
869 
870 
871  if (debug)
872  {
873  // only merging points in debug mode
874 
875  labelList oldToNew;
876  pointField newTgtPoints;
877  bool hasMerged = mergePoints
878  (
879  tgtPoints,
880  SMALL,
881  false,
882  oldToNew,
883  newTgtPoints
884  );
885 
886  if (hasMerged)
887  {
888  Pout<< "Merged from " << tgtPoints.size()
889  << " down to " << newTgtPoints.size() << " points" << endl;
890 
891  tgtPoints.transfer(newTgtPoints);
892  forAll(tgtFaces, i)
893  {
894  inplaceRenumber(oldToNew, tgtFaces[i]);
895  }
896  }
897  }
898 }
899 
900 
901 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::autoPtr::New
static autoPtr< T > New(Args &&... args)
Construct autoPtr of T with forwarding arguments.
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:71
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
InfoInFunction
#define InfoInFunction
Report an information message using Foam::Info.
Definition: messageStream.H:325
Foam::UPstream::commsTypes::nonBlocking
SubField.H
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
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::boundBox::invertedBox
static const boundBox invertedBox
A large inverted boundBox: min/max == +/- ROOTVGREAT.
Definition: boundBox.H:86
meshToMesh.H
Foam::UPstream::parRun
static bool & parRun()
Test if this a parallel run, or allow modify access.
Definition: UPstream.H:434
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:69
Foam::List::append
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:182
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:350
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::inplaceRenumber
void inplaceRenumber(const labelUList &oldToNew, IntListType &input)
Inplace renumber the values (not the indices) of a list.
Definition: ListOpsTemplates.C:61
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
OFstream.H
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Foam::Info
messageStream Info
Information stream (uses stdout - output is on the master only)
DebugInFunction
#define DebugInFunction
Report an information message using Foam::Info.
Definition: messageStream.H:365
Foam::cellList
List< cell > cellList
A List of cells.
Definition: cellListFwd.H:47
AABBTree.H
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::FatalError
error FatalError
processorPolyPatch.H
Foam::add
void add(FieldField< Field1, typename typeOfSum< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
Definition: FieldFieldFunctions.C:939
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::mergePoints
label mergePoints(const PointList &points, const scalar mergeTol, const bool verbose, labelList &pointMap, typename PointList::const_reference origin=PointList::value_type::zero)
Sorts and merges points. All points closer than/equal mergeTol get merged.
Time.H
Foam::autoPtr< Foam::mapDistribute >
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:381
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::allProcs
static rangeType allProcs(const label communicator=worldComm)
Range of process indices for all processes.
Definition: UPstream.H:509
Foam::UPstream::myProcNo
static int myProcNo(const label communicator=worldComm)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:464
Foam::nl
constexpr char nl
Definition: Ostream.H:385
f
labelList f(nPoints)
Foam::faceList
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
mergePoints.H
Merge points. See below.
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
cellBox.H
Foam::DelaunayMeshTools::allPoints
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
cells
const cellShapeList & cells
Definition: gmvOutputHeader.H:3
Foam::treeBoundBoxList
List< treeBoundBox > treeBoundBoxList
List of bounding boxes.
Definition: treeBoundBoxList.H:46
Foam::UPstream::nProcs
static label nProcs(const label communicator=worldComm)
Number of processes in parallel run, and 1 for serial run.
Definition: UPstream.H:446