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