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-2021 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  const auto* procPatch = isA<processorPolyPatch>(pp);
409 
410  // Store info for faces on processor patches
411  const label nbrProci =
412  (procPatch ? procPatch->neighbProcNo() : -1);
413 
414  forAll(pp, i)
415  {
416  label facei = pp.start() + i;
417  label own = tgtMesh.faceOwner()[facei];
418 
419  if (reverseCellMap[own] != -1)
420  {
421  subFaces.append(tgtMesh.faces()[facei]);
422  subFaceOwner.append(reverseCellMap[own]);
423  subFaceNeighbour.append(-1);
424  subNbrProcIDs.append(nbrProci);
425  subProcLocalFaceIDs.append(i);
426  }
427  }
428  }
429 
430  // reverse point map
431  labelList reversePointMap(tgtMesh.nPoints(), -1);
432  DynamicList<point> subPoints(tgtMesh.nPoints());
433  forAll(subFaces, subFacei)
434  {
435  face& f = subFaces[subFacei];
436  forAll(f, fp)
437  {
438  label pointi = f[fp];
439  if (reversePointMap[pointi] == -1)
440  {
441  reversePointMap[pointi] = subPoints.size();
442  subPoints.append(tgtMesh.points()[pointi]);
443  }
444 
445  f[fp] = reversePointMap[pointi];
446  }
447  }
448 
449  // tgt cells into global numbering
450  labelList globalElems(globalI.toGlobal(sendElems));
451 
452  if (debug > 1)
453  {
454  forAll(sendElems, i)
455  {
456  Pout<< "tgtProc:" << Pstream::myProcNo()
457  << " sending tgt cell " << sendElems[i]
458  << "[" << globalElems[i] << "]"
459  << " to srcProc " << domain << endl;
460  }
461  }
462 
463  // pass data
464  if (domain == Pstream::myProcNo())
465  {
466  // allocate my own data
467  points[Pstream::myProcNo()] = subPoints;
468  nInternalFaces[Pstream::myProcNo()] = nInternal;
469  faces[Pstream::myProcNo()] = subFaces;
470  faceOwner[Pstream::myProcNo()] = subFaceOwner;
471  faceNeighbour[Pstream::myProcNo()] = subFaceNeighbour;
472  cellIDs[Pstream::myProcNo()] = globalElems;
473  nbrProcIDs[Pstream::myProcNo()] = subNbrProcIDs;
474  procLocalFaceIDs[Pstream::myProcNo()] = subProcLocalFaceIDs;
475  }
476  else
477  {
478  // send data to other processor domains
479  UOPstream toDomain(domain, pBufs);
480 
481  toDomain
482  << subPoints
483  << nInternal
484  << subFaces
485  << subFaceOwner
486  << subFaceNeighbour
487  << globalElems
488  << subNbrProcIDs
489  << subProcLocalFaceIDs;
490  }
491  }
492  }
493 
494  // Start receiving
495  pBufs.finishedSends();
496 
497  // Consume
498  for (const int domain : Pstream::allProcs())
499  {
500  const labelList& recvElems = map.constructMap()[domain];
501 
502  if (domain != Pstream::myProcNo() && recvElems.size())
503  {
504  UIPstream str(domain, pBufs);
505 
506  str >> points[domain]
507  >> nInternalFaces[domain]
508  >> faces[domain]
509  >> faceOwner[domain]
510  >> faceNeighbour[domain]
511  >> cellIDs[domain]
512  >> nbrProcIDs[domain]
513  >> procLocalFaceIDs[domain];
514  }
515 
516  if (debug)
517  {
518  Pout<< "Target mesh send sizes[" << domain << "]"
519  << ": points="<< points[domain].size()
520  << ", faces=" << faces[domain].size()
521  << ", nInternalFaces=" << nInternalFaces[domain]
522  << ", faceOwn=" << faceOwner[domain].size()
523  << ", faceNbr=" << faceNeighbour[domain].size()
524  << ", cellIDs=" << cellIDs[domain].size() << endl;
525  }
526  }
527 }
528 
529 
530 void Foam::meshToMesh::distributeAndMergeCells
531 (
532  const mapDistribute& map,
533  const polyMesh& tgt,
534  const globalIndex& globalI,
535  pointField& tgtPoints,
536  faceList& tgtFaces,
537  labelList& tgtFaceOwners,
538  labelList& tgtFaceNeighbours,
539  labelList& tgtCellIDs
540 ) const
541 {
542  // Exchange per-processor data
543  List<pointField> allPoints;
544  List<label> allNInternalFaces;
545  List<faceList> allFaces;
546  List<labelList> allFaceOwners;
547  List<labelList> allFaceNeighbours;
548  List<labelList> allTgtCellIDs;
549 
550  // Per target mesh face the neighbouring proc and index in
551  // processor patch (all -1 for normal boundary face)
552  List<labelList> allNbrProcIDs;
553  List<labelList> allProcLocalFaceIDs;
554 
555  distributeCells
556  (
557  map,
558  tgt,
559  globalI,
560  allPoints,
561  allNInternalFaces,
562  allFaces,
563  allFaceOwners,
564  allFaceNeighbours,
565  allTgtCellIDs,
566  allNbrProcIDs,
567  allProcLocalFaceIDs
568  );
569 
570  // Convert lists into format that can be used to generate a valid polyMesh
571  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
572  //
573  // Points and cells are collected into single flat lists:
574  // - i.e. proc0, proc1 ... procN
575  //
576  // Faces need to be sorted after collection to that internal faces are
577  // contiguous, followed by all boundary faces
578  //
579  // Processor patch faces between included cells on neighbouring processors
580  // are converted into internal faces
581  //
582  // Face list structure:
583  // - Per processor:
584  // - internal faces
585  // - processor faces that have been converted into internal faces
586  // - Followed by all boundary faces
587  // - from 'normal' boundary faces
588  // - from singularly-sided processor patch faces
589 
590 
591  // Number of internal+coupled faces
592  labelList allNIntCoupledFaces(allNInternalFaces);
593 
594  // Starting offset for points
595  label nPoints = 0;
596  labelList pointOffset(Pstream::nProcs(), Zero);
597  forAll(allPoints, proci)
598  {
599  pointOffset[proci] = nPoints;
600  nPoints += allPoints[proci].size();
601  }
602 
603  // Starting offset for cells
604  label nCells = 0;
605  labelList cellOffset(Pstream::nProcs(), Zero);
606  forAll(allTgtCellIDs, proci)
607  {
608  cellOffset[proci] = nCells;
609  nCells += allTgtCellIDs[proci].size();
610  }
611 
612  // Count any coupled faces
613  typedef FixedList<label, 3> label3;
614  typedef HashTable<label, label3> procCoupleInfo;
615  procCoupleInfo procFaceToGlobalCell;
616 
617  forAll(allNbrProcIDs, proci)
618  {
619  const labelList& nbrProci = allNbrProcIDs[proci];
620  const labelList& localFacei = allProcLocalFaceIDs[proci];
621 
622  forAll(nbrProci, i)
623  {
624  if (nbrProci[i] != -1 && localFacei[i] != -1)
625  {
626  label3 key;
627  key[0] = min(proci, nbrProci[i]);
628  key[1] = max(proci, nbrProci[i]);
629  key[2] = localFacei[i];
630 
631  const auto fnd = procFaceToGlobalCell.cfind(key);
632 
633  if (!fnd.found())
634  {
635  procFaceToGlobalCell.insert(key, -1);
636  }
637  else
638  {
639  if (debug > 1)
640  {
641  Pout<< "Additional internal face between procs:"
642  << key[0] << " and " << key[1]
643  << " across local face " << key[2] << endl;
644  }
645 
646  allNIntCoupledFaces[proci]++;
647  }
648  }
649  }
650  }
651 
652 
653  // Starting offset for internal faces
654  label nIntFaces = 0;
655  label nFacesTotal = 0;
656  labelList internalFaceOffset(Pstream::nProcs(), Zero);
657  forAll(allNIntCoupledFaces, proci)
658  {
659  label nCoupledFaces =
660  allNIntCoupledFaces[proci] - allNInternalFaces[proci];
661 
662  internalFaceOffset[proci] = nIntFaces;
663  nIntFaces += allNIntCoupledFaces[proci];
664  nFacesTotal += allFaceOwners[proci].size() - nCoupledFaces;
665  }
666 
667  tgtPoints.setSize(nPoints);
668  tgtFaces.setSize(nFacesTotal);
669  tgtFaceOwners.setSize(nFacesTotal);
670  tgtFaceNeighbours.setSize(nFacesTotal);
671  tgtCellIDs.setSize(nCells);
672 
673  // Insert points
674  forAll(allPoints, proci)
675  {
676  const pointField& pts = allPoints[proci];
677  SubList<point>(tgtPoints, pts.size(), pointOffset[proci]) = pts;
678  }
679 
680  // Insert cellIDs
681  forAll(allTgtCellIDs, proci)
682  {
683  const labelList& cellIDs = allTgtCellIDs[proci];
684  SubList<label>(tgtCellIDs, cellIDs.size(), cellOffset[proci]) = cellIDs;
685  }
686 
687 
688  // Insert internal faces (from internal faces)
689  forAll(allFaces, proci)
690  {
691  const faceList& fcs = allFaces[proci];
692  const labelList& faceOs = allFaceOwners[proci];
693  const labelList& faceNs = allFaceNeighbours[proci];
694 
695  SubList<face> slice
696  (
697  tgtFaces,
698  allNInternalFaces[proci],
699  internalFaceOffset[proci]
700  );
701  slice = SubList<face>(fcs, allNInternalFaces[proci]);
702  forAll(slice, i)
703  {
704  add(slice[i], pointOffset[proci]);
705  }
706 
707  SubField<label> ownSlice
708  (
709  tgtFaceOwners,
710  allNInternalFaces[proci],
711  internalFaceOffset[proci]
712  );
713  ownSlice = SubField<label>(faceOs, allNInternalFaces[proci]);
714  add(ownSlice, cellOffset[proci]);
715 
716  SubField<label> nbrSlice
717  (
718  tgtFaceNeighbours,
719  allNInternalFaces[proci],
720  internalFaceOffset[proci]
721  );
722  nbrSlice = SubField<label>(faceNs, allNInternalFaces[proci]);
723  add(nbrSlice, cellOffset[proci]);
724 
725  internalFaceOffset[proci] += allNInternalFaces[proci];
726  }
727 
728 
729  // Insert internal faces (from coupled face-pairs)
730  forAll(allNbrProcIDs, proci)
731  {
732  const labelList& nbrProci = allNbrProcIDs[proci];
733  const labelList& localFacei = allProcLocalFaceIDs[proci];
734  const labelList& faceOs = allFaceOwners[proci];
735  const faceList& fcs = allFaces[proci];
736 
737  forAll(nbrProci, i)
738  {
739  if (nbrProci[i] != -1 && localFacei[i] != -1)
740  {
741  label3 key;
742  key[0] = min(proci, nbrProci[i]);
743  key[1] = max(proci, nbrProci[i]);
744  key[2] = localFacei[i];
745 
746  auto fnd = procFaceToGlobalCell.find(key);
747 
748  if (fnd.found())
749  {
750  label tgtFacei = fnd();
751  if (tgtFacei == -1)
752  {
753  // on first visit store the new cell on this side
754  fnd() = cellOffset[proci] + faceOs[i];
755  }
756  else
757  {
758  // get owner and neighbour in new cell numbering
759  label newOwn = cellOffset[proci] + faceOs[i];
760  label newNbr = fnd();
761  label tgtFacei = internalFaceOffset[proci]++;
762 
763  if (debug > 1)
764  {
765  Pout<< " proc " << proci
766  << "\tinserting face:" << tgtFacei
767  << " connection between owner " << newOwn
768  << " and neighbour " << newNbr
769  << endl;
770  }
771 
772  if (newOwn < newNbr)
773  {
774  // we have correct orientation
775  tgtFaces[tgtFacei] = fcs[i];
776  tgtFaceOwners[tgtFacei] = newOwn;
777  tgtFaceNeighbours[tgtFacei] = newNbr;
778  }
779  else
780  {
781  // reverse orientation
782  tgtFaces[tgtFacei] = fcs[i].reverseFace();
783  tgtFaceOwners[tgtFacei] = newNbr;
784  tgtFaceNeighbours[tgtFacei] = newOwn;
785  }
786 
787  add(tgtFaces[tgtFacei], pointOffset[proci]);
788 
789  // mark with unique value
790  fnd() = -2;
791  }
792  }
793  }
794  }
795  }
796 
797 
798  forAll(allNbrProcIDs, proci)
799  {
800  const labelList& nbrProci = allNbrProcIDs[proci];
801  const labelList& localFacei = allProcLocalFaceIDs[proci];
802  const labelList& faceOs = allFaceOwners[proci];
803  const labelList& faceNs = allFaceNeighbours[proci];
804  const faceList& fcs = allFaces[proci];
805 
806  forAll(nbrProci, i)
807  {
808  // coupled boundary face
809  if (nbrProci[i] != -1 && localFacei[i] != -1)
810  {
811  label3 key;
812  key[0] = min(proci, nbrProci[i]);
813  key[1] = max(proci, nbrProci[i]);
814  key[2] = localFacei[i];
815 
816  label tgtFacei = procFaceToGlobalCell[key];
817 
818  if (tgtFacei == -1)
819  {
820  FatalErrorInFunction
821  << "Unvisited " << key
822  << abort(FatalError);
823  }
824  else if (tgtFacei != -2)
825  {
826  label newOwn = cellOffset[proci] + faceOs[i];
827  label tgtFacei = nIntFaces++;
828 
829  if (debug > 1)
830  {
831  Pout<< " proc " << proci
832  << "\tinserting boundary face:" << tgtFacei
833  << " from coupled face " << key
834  << endl;
835  }
836 
837  tgtFaces[tgtFacei] = fcs[i];
838  add(tgtFaces[tgtFacei], pointOffset[proci]);
839 
840  tgtFaceOwners[tgtFacei] = newOwn;
841  tgtFaceNeighbours[tgtFacei] = -1;
842  }
843  }
844  // normal boundary face
845  else
846  {
847  label own = faceOs[i];
848  label nbr = faceNs[i];
849  if ((own != -1) && (nbr == -1))
850  {
851  label newOwn = cellOffset[proci] + faceOs[i];
852  label tgtFacei = nIntFaces++;
853 
854  tgtFaces[tgtFacei] = fcs[i];
855  add(tgtFaces[tgtFacei], pointOffset[proci]);
856 
857  tgtFaceOwners[tgtFacei] = newOwn;
858  tgtFaceNeighbours[tgtFacei] = -1;
859  }
860  }
861  }
862  }
863 
864 
865  if (debug)
866  {
867  // only merging points in debug mode
868 
869  labelList oldToNew;
870  pointField newTgtPoints;
871  bool hasMerged = mergePoints
872  (
873  tgtPoints,
874  SMALL,
875  false,
876  oldToNew,
877  newTgtPoints
878  );
879 
880  if (hasMerged)
881  {
882  Pout<< "Merged from " << tgtPoints.size()
883  << " down to " << newTgtPoints.size() << " points" << endl;
884 
885  tgtPoints.transfer(newTgtPoints);
886  forAll(tgtFaces, i)
887  {
888  inplaceRenumber(oldToNew, tgtFaces[i]);
889  }
890  }
891  }
892 }
893 
894 
895 // ************************************************************************* //
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:67
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:350
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::glTF::key
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
Definition: foamGltfBase.H:108
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:65
Foam::List::append
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
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 (stdout output on master, null elsewhere)
DebugInFunction
#define DebugInFunction
Report an information message using Foam::Info.
Definition: messageStream.H:388
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: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::allProcs
static rangeType allProcs(const label communicator=worldComm)
Range of process indices for all processes.
Definition: UPstream.H:508
Foam::UPstream::commsTypes::nonBlocking
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
f
labelList f(nPoints)
Foam::faceList
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
Foam::UPstream::parRun
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
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:445