cellsToCellsParallelOps.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) 2012-2023 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 "cellsToCells.H"
27 #include "processorPolyPatch.H"
28 #include "SubField.H"
29 #include "Time.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 inline void offset(label& lst, const label o)
37 {
38  lst += o;
39 }
40 
41 template<class ListType>
42 void offset(ListType& lst, const label o)
43 {
44  forAll(lst, i)
45  {
46  offset(lst[i], o);
47  }
48 }
49 
50 }
51 
52 
53 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
54 
56 (
57  const polyMesh& srcMesh,
58  const polyMesh& tgtMesh
59 ) const
60 {
61  // Calculate and communicate the bound boxes for the source meshes
63  if (srcMesh.nCells() > 0)
64  {
65  // Bounding box for this mesh. Do not reduce.
66  procBb[Pstream::myProcNo()] = boundBox(srcMesh.points(), false);
67 
68  // Slightly increase size to allow for cases where boxes are aligned
69  procBb[Pstream::myProcNo()].inflate(0.01);
70  }
71  Pstream::gatherList(procBb);
72  Pstream::scatterList(procBb);
73 
74  // per processor indices into all segments to send
75  List<DynamicList<label>> resultDyn
76  (
79  );
80 
81  // Send a target cell to a process if it overlaps the source bound box
82  // for that process
83  forAll(tgtMesh.cells(), tgtCelli)
84  {
85  const boundBox cellBb =
86  tgtMesh.cells()[tgtCelli].bb(tgtMesh.points(), tgtMesh.faces());
87 
88  forAll(procBb, proci)
89  {
90  if (procBb[proci].overlaps(cellBb))
91  {
92  resultDyn[proci].append(tgtCelli);
93  }
94  }
95  }
96 
97  // Transfer to permanent storage and return
99  forAll(result, proci)
100  {
101  result[proci].transfer(resultDyn[proci]);
102  }
103 
104  return result;
105 }
106 
107 
109 (
110  const distributionMap& map,
111  const polyMesh& mesh,
112  autoPtr<polyMesh>& localMeshPtr
113 )
114 {
116 
117  // Exchange raw mesh topology/geometry
118  List<pointField> allPoints(Pstream::nProcs());
119  labelList allNInternalFaces(Pstream::nProcs(), 0);
120  List<faceList> allFaces(Pstream::nProcs());
121  List<labelList> allFaceOwners(Pstream::nProcs());
122  List<labelList> allFaceNeighbours(Pstream::nProcs());
123  List<List<remote>> allProcCells(Pstream::nProcs());
124  List<List<remote>> allProcProcessorPatchFaces(Pstream::nProcs());
125  for (label domain = 0; domain < Pstream::nProcs(); domain++)
126  {
127  const labelList& sendElems = map.subMap()[domain];
128 
129  if (sendElems.size())
130  {
131  // reverse cell map
132  labelList reverseCellMap(mesh.nCells(), -1);
133  forAll(sendElems, subCelli)
134  {
135  reverseCellMap[sendElems[subCelli]] = subCelli;
136  }
137 
138  DynamicList<face> subFaces(mesh.nFaces());
139  DynamicList<label> subFaceOwner(mesh.nFaces());
140  DynamicList<label> subFaceNeighbour(mesh.nFaces());
141  DynamicList<remote> subProcProcessorPatchFaces(mesh.nFaces());
142 
143  label nInternal = 0;
144 
145  // internal faces
146  forAll(mesh.faceNeighbour(), facei)
147  {
148  const label own = mesh.faceOwner()[facei];
149  const label nbr = mesh.faceNeighbour()[facei];
150  const label subOwn = reverseCellMap[own];
151  const label subNbr = reverseCellMap[nbr];
152 
153  if (subOwn != -1 && subNbr != -1)
154  {
155  nInternal++;
156 
157  if (subOwn < subNbr)
158  {
159  subFaces.append(mesh.faces()[facei]);
160  subFaceOwner.append(subOwn);
161  subFaceNeighbour.append(subNbr);
162  subProcProcessorPatchFaces.append(remote());
163  }
164  else
165  {
166  subFaces.append(mesh.faces()[facei].reverseFace());
167  subFaceOwner.append(subNbr);
168  subFaceNeighbour.append(subOwn);
169  subProcProcessorPatchFaces.append(remote());
170  }
171  }
172  }
173 
174  // boundary faces for new region
175  forAll(mesh.faceNeighbour(), facei)
176  {
177  const label own = mesh.faceOwner()[facei];
178  const label nbr = mesh.faceNeighbour()[facei];
179  const label subOwn = reverseCellMap[own];
180  const label subNbr = reverseCellMap[nbr];
181 
182  if (subOwn != -1 && subNbr == -1)
183  {
184  subFaces.append(mesh.faces()[facei]);
185  subFaceOwner.append(subOwn);
186  subFaceNeighbour.append(subNbr);
187  subProcProcessorPatchFaces.append(remote());
188  }
189  else if (subOwn == -1 && subNbr != -1)
190  {
191  subFaces.append(mesh.faces()[facei].reverseFace());
192  subFaceOwner.append(subNbr);
193  subFaceNeighbour.append(subOwn);
194  subProcProcessorPatchFaces.append(remote());
195  }
196  }
197 
198  // boundary faces of existing region
199  forAll(mesh.boundaryMesh(), patchi)
200  {
201  const polyPatch& pp = mesh.boundaryMesh()[patchi];
202 
203  const label nbrProci =
204  isType<processorPolyPatch>(pp)
205  ? refCast<const processorPolyPatch>(pp).neighbProcNo()
206  : -1;
207 
208  forAll(pp, i)
209  {
210  const label facei = pp.start() + i;
211  const label own = mesh.faceOwner()[facei];
212 
213  if (reverseCellMap[own] != -1)
214  {
215  subFaces.append(mesh.faces()[facei]);
216  subFaceOwner.append(reverseCellMap[own]);
217  subFaceNeighbour.append(-1);
218  subProcProcessorPatchFaces.append
219  (
220  nbrProci != -1 ? remote(nbrProci, i) : remote()
221  );
222  }
223  }
224  }
225 
226  // reverse point map
227  labelList reversePointMap(mesh.nPoints(), -1);
228  DynamicList<point> subPoints(mesh.nPoints());
229  forAll(subFaces, subFacei)
230  {
231  face& f = subFaces[subFacei];
232  forAll(f, fp)
233  {
234  label pointi = f[fp];
235  if (reversePointMap[pointi] == -1)
236  {
237  reversePointMap[pointi] = subPoints.size();
238  subPoints.append(mesh.points()[pointi]);
239  }
240 
241  f[fp] = reversePointMap[pointi];
242  }
243  }
244 
245  // cell indices
246  List<remote> subProcCells(sendElems.size());
247  forAll(sendElems, i)
248  {
249  subProcCells[i] = remote(Pstream::myProcNo(), sendElems[i]);
250  }
251 
252  // pass data
253  if (domain == Pstream::myProcNo())
254  {
255  // allocate my own data
256  allPoints[Pstream::myProcNo()] = subPoints;
257  allNInternalFaces[Pstream::myProcNo()] = nInternal;
258  allFaces[Pstream::myProcNo()] = subFaces;
259  allFaceOwners[Pstream::myProcNo()] = subFaceOwner;
260  allFaceNeighbours[Pstream::myProcNo()] = subFaceNeighbour;
261  allProcCells[Pstream::myProcNo()] = subProcCells;
262  allProcProcessorPatchFaces[Pstream::myProcNo()] =
263  subProcProcessorPatchFaces;
264  }
265  else
266  {
267  // send data to other processor domains
268  UOPstream toDomain(domain, pBufs);
269 
270  toDomain
271  << subPoints
272  << nInternal
273  << subFaces
274  << subFaceOwner
275  << subFaceNeighbour
276  << subProcCells
277  << subProcProcessorPatchFaces;
278  }
279  }
280  }
281 
282  // Start receiving
283  pBufs.finishedSends();
284 
285  // Consume
286  for (label domain = 0; domain < Pstream::nProcs(); domain++)
287  {
288  const labelList& recvElems = map.constructMap()[domain];
289 
290  if (domain != Pstream::myProcNo() && recvElems.size())
291  {
292  UIPstream str(domain, pBufs);
293 
294  str >> allPoints[domain]
295  >> allNInternalFaces[domain]
296  >> allFaces[domain]
297  >> allFaceOwners[domain]
298  >> allFaceNeighbours[domain]
299  >> allProcCells[domain]
300  >> allProcProcessorPatchFaces[domain];
301  }
302  }
303 
304  // Convert lists into format that can be used to generate a valid polyMesh
305  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
306  //
307  // Points and cells are collected into single flat lists:
308  // - i.e. proc0, proc1 ... procN
309  //
310  // Faces need to be sorted after collection to that internal faces are
311  // contiguous, followed by all boundary faces
312  //
313  // Processor patch faces between included cells on neighbouring processors
314  // are converted into internal faces
315  //
316  // Face list structure:
317  // - Per processor:
318  // - internal faces
319  // - processor faces that have been converted into internal faces
320  // - Followed by all boundary faces
321  // - from 'normal' boundary faces
322  // - from singularly-sided processor patch faces
323 
324  // Number of internal+coupled faces
325  labelList allNIntCoupledFaces(allNInternalFaces);
326 
327  // Starting offset for points
328  label nPoints = 0;
329  labelList pointOffset(Pstream::nProcs(), 0);
330  forAll(allPoints, proci)
331  {
332  pointOffset[proci] = nPoints;
333  nPoints += allPoints[proci].size();
334  }
335 
336  // Count any coupled faces
337  typedef FixedList<label, 3> label3;
338  typedef HashTable<label, label3, label3::Hash<>> procCoupleInfo;
339  procCoupleInfo procFaceToGlobalCell;
340  forAll(allProcProcessorPatchFaces, proci)
341  {
342  const List<remote>& procProcessorPatchFaces =
343  allProcProcessorPatchFaces[proci];
344 
345  forAll(procProcessorPatchFaces, i)
346  {
347  if (procProcessorPatchFaces[i] != remote())
348  {
349  const label3 key
350  ({
351  min(proci, procProcessorPatchFaces[i].proci),
352  max(proci, procProcessorPatchFaces[i].proci),
353  procProcessorPatchFaces[i].elementi
354  });
355 
356  procCoupleInfo::const_iterator fnd =
357  procFaceToGlobalCell.find(key);
358 
359  if (fnd == procFaceToGlobalCell.end())
360  {
361  procFaceToGlobalCell.insert(key, -1);
362  }
363  else
364  {
365  if (debug)
366  {
367  Pout<< "Additional internal face between procs:"
368  << key[0] << " and " << key[1]
369  << " across local face " << key[2] << endl;
370  }
371 
372  allNIntCoupledFaces[proci]++;
373  }
374  }
375  }
376  }
377 
378  // Starting offset for internal faces
379  label nIntFaces = 0;
380  label nFacesTotal = 0;
381  labelList internalFaceOffset(Pstream::nProcs(), 0);
382  forAll(allNIntCoupledFaces, proci)
383  {
384  label nCoupledFaces =
385  allNIntCoupledFaces[proci] - allNInternalFaces[proci];
386 
387  internalFaceOffset[proci] = nIntFaces;
388  nIntFaces += allNIntCoupledFaces[proci];
389  nFacesTotal += allFaceOwners[proci].size() - nCoupledFaces;
390  }
391 
392  // Starting offset for cells
393  label nCells = 0;
394  labelList cellOffset(Pstream::nProcs(), 0);
395  forAll(allProcCells, proci)
396  {
397  cellOffset[proci] = nCells;
398  nCells += allProcCells[proci].size();
399  }
400 
401  // Allocate
402  List<remote> localProcCells(nCells);
403  pointField localPoints(nPoints);
404  faceList localFaces(nFacesTotal);
405  labelList localFaceOwners(nFacesTotal);
406  labelList localFaceNeighbours(nFacesTotal);
407 
408  // Insert proc-cells
409  forAll(allProcCells, proci)
410  {
411  const List<remote>& procCells = allProcCells[proci];
412  SubList<remote>(localProcCells, procCells.size(), cellOffset[proci]) =
413  procCells;
414  }
415 
416  // Insert points
417  forAll(allPoints, proci)
418  {
419  const pointField& pts = allPoints[proci];
420  SubList<point>(localPoints, pts.size(), pointOffset[proci]) = pts;
421  }
422 
423  // Insert internal faces (from internal faces)
424  forAll(allFaces, proci)
425  {
426  const faceList& fcs = allFaces[proci];
427  const labelList& faceOs = allFaceOwners[proci];
428  const labelList& faceNs = allFaceNeighbours[proci];
429 
430  SubList<face> slice
431  (
432  localFaces,
433  allNInternalFaces[proci],
434  internalFaceOffset[proci]
435  );
436  slice = SubList<face>(fcs, allNInternalFaces[proci]);
437  offset(slice, pointOffset[proci]);
438 
439  SubField<label> ownSlice
440  (
441  localFaceOwners,
442  allNInternalFaces[proci],
443  internalFaceOffset[proci]
444  );
445  ownSlice = SubField<label>(faceOs, allNInternalFaces[proci]);
446  offset(ownSlice, cellOffset[proci]);
447 
448  SubField<label> nbrSlice
449  (
450  localFaceNeighbours,
451  allNInternalFaces[proci],
452  internalFaceOffset[proci]
453  );
454  nbrSlice = SubField<label>(faceNs, allNInternalFaces[proci]);
455  offset(nbrSlice, cellOffset[proci]);
456 
457  internalFaceOffset[proci] += allNInternalFaces[proci];
458  }
459 
460  // Insert internal faces (from coupled face-pairs)
461  forAll(allProcProcessorPatchFaces, proci)
462  {
463  const List<remote>& procProcessorPatchFaces =
464  allProcProcessorPatchFaces[proci];
465  const labelList& faceOs = allFaceOwners[proci];
466  const faceList& fcs = allFaces[proci];
467 
468  forAll(procProcessorPatchFaces, i)
469  {
470  if (procProcessorPatchFaces[i] != remote())
471  {
472  const label3 key
473  ({
474  min(proci, procProcessorPatchFaces[i].proci),
475  max(proci, procProcessorPatchFaces[i].proci),
476  procProcessorPatchFaces[i].elementi
477  });
478 
479  procCoupleInfo::iterator fnd = procFaceToGlobalCell.find(key);
480 
481  if (fnd != procFaceToGlobalCell.end())
482  {
483  if (fnd() == -1)
484  {
485  // on first visit store the new cell on this side
486  fnd() = cellOffset[proci] + faceOs[i];
487  }
488  else
489  {
490  // get owner and neighbour in new cell numbering
491  const label newOwn = cellOffset[proci] + faceOs[i];
492  const label newNbr = fnd();
493  const label localFacei = internalFaceOffset[proci]++;
494 
495  if (debug)
496  {
497  Pout<< " proc " << proci
498  << "\tinserting face:" << localFacei
499  << " connection between owner " << newOwn
500  << " and neighbour " << newNbr
501  << endl;
502  }
503 
504  if (newOwn < newNbr)
505  {
506  // we have correct orientation
507  localFaces[localFacei] = fcs[i];
508  localFaceOwners[localFacei] = newOwn;
509  localFaceNeighbours[localFacei] = newNbr;
510  }
511  else
512  {
513  // reverse orientation
514  localFaces[localFacei] = fcs[i].reverseFace();
515  localFaceOwners[localFacei] = newNbr;
516  localFaceNeighbours[localFacei] = newOwn;
517  }
518 
519  offset(localFaces[localFacei], pointOffset[proci]);
520 
521  // mark with unique value
522  fnd() = -2;
523  }
524  }
525  }
526  }
527  }
528 
529  forAll(allProcProcessorPatchFaces, proci)
530  {
531  const List<remote>& procProcessorPatchFaces =
532  allProcProcessorPatchFaces[proci];
533  const labelList& faceOs = allFaceOwners[proci];
534  const labelList& faceNs = allFaceNeighbours[proci];
535  const faceList& fcs = allFaces[proci];
536 
537  forAll(procProcessorPatchFaces, i)
538  {
539  // coupled boundary face
540  if (procProcessorPatchFaces[i] != remote())
541  {
542  const label3 key
543  ({
544  min(proci, procProcessorPatchFaces[i].proci),
545  max(proci, procProcessorPatchFaces[i].proci),
546  procProcessorPatchFaces[i].elementi
547  });
548 
549  if (procFaceToGlobalCell[key] == -1)
550  {
552  << "Unvisited " << key
553  << abort(FatalError);
554  }
555  else if (procFaceToGlobalCell[key] != -2)
556  {
557  const label newOwn = cellOffset[proci] + faceOs[i];
558  const label localFacei = nIntFaces++;
559 
560  if (debug)
561  {
562  Pout<< " proc " << proci
563  << "\tinserting boundary face:" << localFacei
564  << " from coupled face " << key
565  << endl;
566  }
567 
568  localFaces[localFacei] = fcs[i];
569  offset(localFaces[localFacei], pointOffset[proci]);
570 
571  localFaceOwners[localFacei] = newOwn;
572  localFaceNeighbours[localFacei] = -1;
573  }
574  }
575 
576  // normal boundary face
577  else
578  {
579  const label own = faceOs[i];
580  const label nbr = faceNs[i];
581 
582  if ((own != -1) && (nbr == -1))
583  {
584  const label newOwn = cellOffset[proci] + faceOs[i];
585  const label localFacei = nIntFaces++;
586 
587  localFaces[localFacei] = fcs[i];
588  offset(localFaces[localFacei], pointOffset[proci]);
589 
590  localFaceOwners[localFacei] = newOwn;
591  localFaceNeighbours[localFacei] = -1;
592  }
593  }
594  }
595  }
596 
597  // Create the local mesh
598  localMeshPtr.reset
599  (
600  new polyMesh
601  (
602  IOobject
603  (
604  "local" + mesh.name().capitalise(),
605  mesh.time().name(),
606  mesh.time(),
608  ),
609  move(localPoints),
610  move(localFaces),
611  move(localFaceOwners),
612  move(localFaceNeighbours),
613  false
614  )
615  );
616 
617  // Add a dummy patch to the target mesh
619  patches[0] = new polyPatch
620  (
621  "defaultFaces",
622  localMeshPtr().nFaces() - localMeshPtr().nInternalFaces(),
623  localMeshPtr().nInternalFaces(),
624  0,
625  localMeshPtr().boundaryMesh(),
626  word::null
627  );
628  localMeshPtr().addPatches(patches);
629 
630  // Force calculation of tet-base points used for point-in-cell
631  (void) localMeshPtr().tetBasePtIs();
632 
633  return localProcCells;
634 }
635 
636 
638 {
639  // Determine which local target cells are actually used
640  boolList oldLocalTgtCellIsUsed(localTgtProcCellsPtr_().size(), false);
641  forAll(srcLocalTgtCells_, srcCelli)
642  {
643  forAll(srcLocalTgtCells_[srcCelli], i)
644  {
645  oldLocalTgtCellIsUsed[srcLocalTgtCells_[srcCelli][i]] = true;
646  }
647  }
648 
649  // Trim the target map
650  labelList oldToNewLocalTgtCell, newToOldLocalTgtCell;
652  (
653  oldLocalTgtCellIsUsed,
654  tgtMapPtr_(),
655  oldToNewLocalTgtCell,
656  newToOldLocalTgtCell
657  );
658 
659  if (debug)
660  {
661  Pout<< "Trim from " << oldToNewLocalTgtCell.size() << " to "
662  << newToOldLocalTgtCell.size() << " cells" << endl;
663  }
664 
665  // Renumber the source addressing
666  forAll(srcLocalTgtCells_, srcCelli)
667  {
668  forAll(srcLocalTgtCells_[srcCelli], i)
669  {
670  srcLocalTgtCells_[srcCelli][i] =
671  oldToNewLocalTgtCell[srcLocalTgtCells_[srcCelli][i]];
672  }
673  }
674 
675  // Trim the target addressing
676  tgtLocalSrcCells_ =
677  labelListList(tgtLocalSrcCells_, newToOldLocalTgtCell);
678  localTgtProcCellsPtr_() =
679  List<remote>(localTgtProcCellsPtr_(), newToOldLocalTgtCell);
680 
681  // Trim the target weights
682  tgtWeights_ = scalarListList(tgtWeights_, newToOldLocalTgtCell);
683 
684  // Trim the stored local target mesh
685  const polyMesh& oldLocalTgtMesh = localTgtMeshPtr_();
686  const labelList& oldLocalTgtFaceOwner =
687  oldLocalTgtMesh.faceOwner();
688  labelList oldLocalTgtFaceNeighbour(oldLocalTgtFaceOwner.size(), -1);
689  SubList<label>(oldLocalTgtFaceNeighbour, oldLocalTgtMesh.nInternalFaces()) =
690  oldLocalTgtMesh.faceNeighbour();
691 
692  // ...
693  labelList newToOldLocalTgtFace(identityMap(oldLocalTgtMesh.nFaces()));
694  labelList oldToNewLocalTgtFace;
695  {
696  label i0 = 0;
697  label i1 = newToOldLocalTgtFace.size();
698  label iEnd = newToOldLocalTgtFace.size();
699 
700  while (i0 < i1)
701  {
702  label& oldLocalTgtFacei0 = newToOldLocalTgtFace[i0];
703  label& oldLocalTgtFacei1 = newToOldLocalTgtFace[i1 - 1];
704  label& oldLocalTgtFaceiEnd = newToOldLocalTgtFace[iEnd - 1];
705 
706  const label newLocalTgtOwni0 =
707  oldLocalTgtFaceOwner[oldLocalTgtFacei0] != -1
708  ? oldToNewLocalTgtCell
709  [oldLocalTgtFaceOwner[oldLocalTgtFacei0]]
710  : -1;
711  const label newLocalTgtOwni1 =
712  oldLocalTgtFaceOwner[oldLocalTgtFacei1] != -1
713  ? oldToNewLocalTgtCell
714  [oldLocalTgtFaceOwner[oldLocalTgtFacei1]]
715  : -1;
716 
717  const label newLocalTgtNbri0 =
718  oldLocalTgtFaceNeighbour[oldLocalTgtFacei0] != -1
719  ? oldToNewLocalTgtCell
720  [oldLocalTgtFaceNeighbour[oldLocalTgtFacei0]]
721  : -1;
722  const label newLocalTgtNbri1 =
723  oldLocalTgtFaceNeighbour[oldLocalTgtFacei1] != -1
724  ? oldToNewLocalTgtCell
725  [oldLocalTgtFaceNeighbour[oldLocalTgtFacei1]]
726  : -1;
727 
728  const bool used0 =
729  newLocalTgtOwni0 != -1 || newLocalTgtNbri0 != -1;
730  const bool used1 =
731  newLocalTgtOwni1 != -1 || newLocalTgtNbri1 != -1;
732 
733  const bool internal0 =
734  newLocalTgtOwni0 != -1 && newLocalTgtNbri0 != -1;
735  const bool internal1 =
736  newLocalTgtOwni1 != -1 && newLocalTgtNbri1 != -1;
737 
738  // If face 0 is not used, move it to the end to remove it
739  if (!used0)
740  {
741  Swap(oldLocalTgtFacei0, oldLocalTgtFaceiEnd);
742  if (i1 == iEnd) i1 --;
743  iEnd --;
744  }
745 
746  // If face 1 is not used, move it to the end to remove it
747  else if (!used1)
748  {
749  Swap(oldLocalTgtFacei1, oldLocalTgtFaceiEnd);
750  if (i1 == iEnd) i1 --;
751  iEnd --;
752  }
753 
754  // Both are internal faces. Face 0 is fine, but face 1 might be out
755  // of order. So move to the next face 0.
756  else if (internal0 && internal1)
757  {
758  i0 ++;
759  }
760 
761  // Both are boundary faces. Face 0 might be out of order, but face
762  // 1 is fine. So move to the next face 1.
763  else if (!internal0 && !internal1)
764  {
765  i1 --;
766  }
767 
768  // Face 0 is an internal face and face 1 is a boundary face. Both
769  // are fine. So move to the next two faces.
770  else if (internal0 && !internal1)
771  {
772  i0 ++;
773  i1 --;
774  }
775 
776  // Face 0 is a boundary face and face 1 is an internal face. Both
777  // are out of order. So swap and then move to the next two faces.
778  else if (!internal0 && internal1)
779  {
780  Swap(oldLocalTgtFacei0, oldLocalTgtFacei1);
781  i0 ++;
782  i1 --;
783  }
784  }
785 
786  newToOldLocalTgtFace.resize(iEnd);
787 
788  oldToNewLocalTgtFace =
789  invert(oldLocalTgtMesh.nFaces(), newToOldLocalTgtFace);
790  }
791 
792  if (debug)
793  {
794  Pout<< "Trim from " << oldToNewLocalTgtFace.size() << " to "
795  << newToOldLocalTgtFace.size() << " faces" << endl;
796  }
797 
798  // Create trimmed mesh primitives
799  pointField newLocalTgtPoints(oldLocalTgtMesh.points());
800  faceList newLocalTgtFaces(oldLocalTgtMesh.faces(), newToOldLocalTgtFace);
801  labelList newLocalTgtFaceOwner
802  (
803  oldLocalTgtFaceOwner,
804  static_cast<const labelUList&>(newToOldLocalTgtFace)
805  );
806  inplaceRenumber(oldToNewLocalTgtCell, newLocalTgtFaceOwner);
807  labelList newLocalTgtFaceNeighbour
808  (
809  oldLocalTgtFaceNeighbour,
810  static_cast<const labelUList&>(newToOldLocalTgtFace)
811  );
812  inplaceRenumber(oldToNewLocalTgtCell, newLocalTgtFaceNeighbour);
813 
814  // Check the trimmed mesh structure and flip any faces that have a
815  // neighbour but not an owner
816  {
817  label newLocalTgtNInternalFaces = 0;
818  bool internal0 = true;
819 
820  forAll(newLocalTgtFaces, newLocalTgtFacei)
821  {
822  face& newLocalTgtF = newLocalTgtFaces[newLocalTgtFacei];
823  label& newLocalTgtOwni = newLocalTgtFaceOwner[newLocalTgtFacei];
824  label& newLocalTgtNbri = newLocalTgtFaceNeighbour[newLocalTgtFacei];
825 
826  // Check that internal and boundary faces are in order
827  const bool internal =
828  newLocalTgtOwni != -1 && newLocalTgtNbri != -1;
829  if (internal0 && !internal)
830  {
831  newLocalTgtNInternalFaces = newLocalTgtFacei;
832  internal0 = false;
833  }
834  if (!internal0 && internal)
835  {
837  << "Trimmed mesh has boundary faces before internal faces"
838  << exit(FatalError);
839  }
840 
841  // Flip any face that has a neighbour but not an owner
842  const bool flip =
843  newLocalTgtOwni == -1 && newLocalTgtNbri != -1;
844  if (flip)
845  {
846  newLocalTgtF = newLocalTgtF.reverseFace();
847  Swap(newLocalTgtOwni, newLocalTgtNbri);
848  }
849  }
850 
851  newLocalTgtFaceNeighbour.resize(newLocalTgtNInternalFaces);
852  }
853 
854  // Create the local mesh
855  localTgtMeshPtr_.reset
856  (
857  new polyMesh
858  (
859  IOobject
860  (
861  "trimmed" + oldLocalTgtMesh.name().capitalise(),
862  oldLocalTgtMesh.time().name(),
863  oldLocalTgtMesh.time(),
865  ),
866  move(newLocalTgtPoints),
867  move(newLocalTgtFaces),
868  move(newLocalTgtFaceOwner),
869  move(newLocalTgtFaceNeighbour),
870  false
871  )
872  );
873 
874  // Add a dummy patch to the target mesh
876  patches[0] = new polyPatch
877  (
878  "defaultFaces",
879  localTgtMeshPtr_().nFaces() - localTgtMeshPtr_().nInternalFaces(),
880  localTgtMeshPtr_().nInternalFaces(),
881  0,
882  localTgtMeshPtr_().boundaryMesh(),
883  word::null
884  );
885  localTgtMeshPtr_().addPatches(patches);
886 
887  // Force calculation of tet-base points used for point-in-cell
888  (void) localTgtMeshPtr_().tetBasePtIs();
889 }
890 
891 
892 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
An STL-conforming hash table.
Definition: HashTable.H:127
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
const word & name() const
Return name.
Definition: IOobject.H:310
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
void resize(const label)
Alias for setSize(const label)
Definition: ListI.H:138
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Buffers for inter-processor communications streams (UOPstream, UIPstream).
void finishedSends(const bool block=true)
Mark all sends as having been done. This will start receives.
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
Pre-declare related SubField type.
Definition: SubField.H:63
A List obtained as a section of another List.
Definition: SubList.H:56
Input inter-processor communications stream operating on external buffer.
Definition: UIPstream.H:57
Output inter-processor communications stream operating on external buffer.
Definition: UOPstream.H:58
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: autoPtr.H:51
void reset(T *=nullptr)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:59
labelListList tgtMeshSendCells(const polyMesh &srcMesh, const polyMesh &tgtMesh) const
Determine which target cells need to be sent to the source.
void trimLocalTgt()
Trim the local target addressing and mesh so that communication.
static List< remote > distributeMesh(const distributionMap &map, const polyMesh &mesh, autoPtr< polyMesh > &localMeshPtr)
Distribute a mesh given its distribution map.
const word & name() const
Return const reference to name.
const labelListList & constructMap() const
From subsetted data to new reconstructed data.
const labelListList & subMap() const
From subsetted data back to original data.
Class containing processor-to-processor mapping information.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:76
face reverseFace() const
Return face with reverse direction.
Definition: face.C:256
const Time & time() const
Return time.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1326
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1339
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:404
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1345
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1313
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:70
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:280
label nCells() const
label nPoints() const
label nInternalFaces() const
label nFaces() const
const cellList & cells() const
Struct for keeping processor, element (cell, face, point) index.
Definition: remote.H:57
word capitalise() const
Return the word with the first letter capitalised.
Definition: wordI.H:131
static const word null
An empty word.
Definition: word.H:77
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
label patchi
label nPoints
const fvPatchList & patches
void trimDistributionMap(const boolList &oldIsUsed, distributionMap &map, labelList &oldToNew, labelList &newToOld)
Trim unused elements from a distribution map. Return the map and maps.
Namespace for OpenFOAM.
List< scalarList > scalarListList
Definition: scalarList.H:51
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
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:59
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
Definition: ListOps.C:37
errorManip< error > abort(error &err)
Definition: errorManip.H:131
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
error FatalError
void Swap(T &a, T &b)
Definition: Swap.H:43
labelList identityMap(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
void offset(label &lst, const label o)
labelList f(nPoints)