fvMeshDistribute.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) 2011-2018 OpenFOAM Foundation
9  Copyright (C) 2015-2022 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 "fvMeshDistribute.H"
30 #include "fvMeshAdder.H"
31 #include "fvMeshSubset.H"
32 #include "faceCoupleInfo.H"
33 #include "processorFvPatchField.H"
34 #include "processorFvsPatchField.H"
35 #include "processorCyclicPolyPatch.H"
36 #include "processorCyclicFvPatchField.H"
37 #include "polyTopoChange.H"
38 #include "removeCells.H"
39 #include "polyModifyFace.H"
40 #include "polyRemovePoint.H"
41 #include "mapDistributePolyMesh.H"
42 #include "surfaceFields.H"
43 #include "syncTools.H"
44 #include "CompactListList.H"
45 #include "fvMeshTools.H"
46 #include "labelPairHashes.H"
47 #include "ListOps.H"
48 #include "globalIndex.H"
49 #include "cyclicACMIPolyPatch.H"
50 #include "mappedPatchBase.H"
51 
52 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
53 
54 namespace Foam
55 {
56  defineTypeNameAndDebug(fvMeshDistribute, 0);
57 
58  //- Less function class that can be used for sorting processor patches
59  class lessProcPatches
60  {
61  const labelList& nbrProc_;
62  const labelList& referPatchID_;
63 
64  public:
65 
66  lessProcPatches(const labelList& nbrProc, const labelList& referPatchID)
67  :
68  nbrProc_(nbrProc),
69  referPatchID_(referPatchID)
70  {}
71 
72  bool operator()(const label a, const label b)
73  {
74  if (nbrProc_[a] < nbrProc_[b])
75  {
76  return true;
77  }
78  else if (nbrProc_[a] > nbrProc_[b])
79  {
80  return false;
81  }
82  else
83  {
84  // Equal neighbour processor
85  return referPatchID_[a] < referPatchID_[b];
86  }
87  }
88  };
89 }
90 
91 
92 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
93 
94 void Foam::fvMeshDistribute::inplaceRenumberWithFlip
95 (
96  const labelUList& oldToNew,
97  const bool oldToNewHasFlip,
98  const bool lstHasFlip,
99  labelUList& lst
100 )
101 {
102  if (!lstHasFlip && !oldToNewHasFlip)
103  {
104  Foam::inplaceRenumber(oldToNew, lst);
105  }
106  else
107  {
108  // Either input data or map encodes sign so result encodes sign
109 
110  forAll(lst, elemI)
111  {
112  // Extract old value and sign
113  label val = lst[elemI];
114  label sign = 1;
115  if (lstHasFlip)
116  {
117  if (val > 0)
118  {
119  val = val-1;
120  }
121  else if (val < 0)
122  {
123  val = -val-1;
124  sign = -1;
125  }
126  else
127  {
128  FatalErrorInFunction
129  << "Problem : zero value " << val
130  << " at index " << elemI << " out of " << lst.size()
131  << " list with flip bit" << exit(FatalError);
132  }
133  }
134 
135 
136  // Lookup new value and possibly change sign
137  label newVal = oldToNew[val];
138 
139  if (oldToNewHasFlip)
140  {
141  if (newVal > 0)
142  {
143  newVal = newVal-1;
144  }
145  else if (newVal < 0)
146  {
147  newVal = -newVal-1;
148  sign = -sign;
149  }
150  else
151  {
152  FatalErrorInFunction
153  << "Problem : zero value " << newVal
154  << " at index " << elemI << " out of "
155  << oldToNew.size()
156  << " list with flip bit" << exit(FatalError);
157  }
158  }
159 
160 
161  // Encode new value and sign
162  lst[elemI] = sign*(newVal+1);
163  }
164  }
165 }
166 
167 
168 Foam::labelList Foam::fvMeshDistribute::select
169 (
170  const bool selectEqual,
171  const labelList& values,
172  const label value
173 )
174 {
175  label n = 0;
176 
177  forAll(values, i)
178  {
179  if (selectEqual == (values[i] == value))
180  {
181  n++;
182  }
183  }
184 
185  labelList indices(n);
186  n = 0;
187 
188  forAll(values, i)
189  {
190  if (selectEqual == (values[i] == value))
191  {
192  indices[n++] = i;
193  }
194  }
195  return indices;
196 }
197 
198 
199 Foam::wordList Foam::fvMeshDistribute::mergeWordList(const wordList& procNames)
200 {
201  List<wordList> allNames(Pstream::nProcs());
202  allNames[Pstream::myProcNo()] = procNames;
203  Pstream::gatherList(allNames);
204 
205  // Assume there are few zone names to merge. Use HashSet otherwise (but
206  // maintain ordering)
207  DynamicList<word> mergedNames;
208  if (Pstream::master())
209  {
210  mergedNames = procNames;
211  for (const wordList& names : allNames)
212  {
213  for (const word& name : names)
214  {
215  mergedNames.appendUniq(name);
216  }
217  }
218  }
219  Pstream::broadcast(mergedNames);
220 
221  return mergedNames;
222 }
223 
224 
225 void Foam::fvMeshDistribute::printMeshInfo(const fvMesh& mesh)
226 {
227  Pout<< "Primitives:" << nl
228  << " points :" << mesh.nPoints() << nl
229  << " bb :" << boundBox(mesh.points(), false) << nl
230  << " internalFaces:" << mesh.nInternalFaces() << nl
231  << " faces :" << mesh.nFaces() << nl
232  << " cells :" << mesh.nCells() << nl;
233 
234  const fvBoundaryMesh& patches = mesh.boundary();
235 
236  Pout<< "Patches:" << endl;
237  forAll(patches, patchi)
238  {
239  const polyPatch& pp = patches[patchi].patch();
240 
241  Pout<< " " << patchi << " name:" << pp.name()
242  << " size:" << pp.size()
243  << " start:" << pp.start()
244  << " type:" << pp.type()
245  << endl;
246  }
247 
248  if (mesh.pointZones().size())
249  {
250  Pout<< "PointZones:" << endl;
251  forAll(mesh.pointZones(), zoneI)
252  {
253  const pointZone& pz = mesh.pointZones()[zoneI];
254  Pout<< " " << zoneI << " name:" << pz.name()
255  << " size:" << pz.size()
256  << endl;
257  }
258  }
259  if (mesh.faceZones().size())
260  {
261  Pout<< "FaceZones:" << endl;
262  forAll(mesh.faceZones(), zoneI)
263  {
264  const faceZone& fz = mesh.faceZones()[zoneI];
265  Pout<< " " << zoneI << " name:" << fz.name()
266  << " size:" << fz.size()
267  << endl;
268  }
269  }
270  if (mesh.cellZones().size())
271  {
272  Pout<< "CellZones:" << endl;
273  forAll(mesh.cellZones(), zoneI)
274  {
275  const cellZone& cz = mesh.cellZones()[zoneI];
276  Pout<< " " << zoneI << " name:" << cz.name()
277  << " size:" << cz.size()
278  << endl;
279  }
280  }
281 }
282 
283 
284 void Foam::fvMeshDistribute::printCoupleInfo
285 (
286  const primitiveMesh& mesh,
287  const labelList& sourceFace,
288  const labelList& sourceProc,
289  const labelList& sourcePatch,
290  const labelList& sourceNewNbrProc
291 )
292 {
293  Pout<< nl
294  << "Current coupling info:"
295  << endl;
296 
297  forAll(sourceFace, bFacei)
298  {
299  label meshFacei = mesh.nInternalFaces() + bFacei;
300 
301  Pout<< " meshFace:" << meshFacei
302  << " fc:" << mesh.faceCentres()[meshFacei]
303  << " connects to proc:" << sourceProc[bFacei]
304  << "/face:" << sourceFace[bFacei]
305  << " which will move to proc:" << sourceNewNbrProc[bFacei]
306  << endl;
307  }
308 }
309 
310 
311 Foam::label Foam::fvMeshDistribute::findNonEmptyPatch() const
312 {
313  // Finds (non-empty) patch that exposed internal and proc faces can be
314  // put into.
315  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
316 
317 
318  // Mark 'special' patches : -coupled, -duplicate faces. These probably
319  // should not be used to (temporarily) store processor faces ...
320 
321  bitSet isCoupledPatch(patches.size());
322  forAll(patches, patchi)
323  {
324  const polyPatch& pp = patches[patchi];
325  const auto* cpp = isA<cyclicACMIPolyPatch>(pp);
326 
327  if (cpp)
328  {
329  isCoupledPatch.set(patchi);
330  const label dupPatchID = cpp->nonOverlapPatchID();
331  if (dupPatchID != -1)
332  {
333  isCoupledPatch.set(dupPatchID);
334  }
335  }
336  else if (pp.coupled())
337  {
338  isCoupledPatch.set(patchi);
339  }
340  }
341 
342  label nonEmptyPatchi = -1;
343 
344  forAllReverse(patches, patchi)
345  {
346  const polyPatch& pp = patches[patchi];
347 
348  if
349  (
350  !isA<emptyPolyPatch>(pp)
351  && !isCoupledPatch(patchi)
352  && !isA<mappedPatchBase>(pp)
353  )
354  {
355  nonEmptyPatchi = patchi;
356  break;
357  }
358  }
359 
360  if (nonEmptyPatchi == -1)
361  {
362  FatalErrorInFunction
363  << "Cannot find a patch which is not of type empty, mapped or"
364  << " coupled in patches " << patches.names() << endl
365  << "There has to be at least one such patch for"
366  << " distribution to work" << abort(FatalError);
367  }
368 
369  if (debug)
370  {
371  Pout<< "findNonEmptyPatch : using patch " << nonEmptyPatchi
372  << " name:" << patches[nonEmptyPatchi].name()
373  << " type:" << patches[nonEmptyPatchi].type()
374  << " to put exposed faces into." << endl;
375  }
376 
377 
378  // Do additional test for processor patches intermingled with non-proc
379  // patches.
380  label procPatchi = -1;
381 
382  forAll(patches, patchi)
383  {
384  if (isA<processorPolyPatch>(patches[patchi]))
385  {
386  procPatchi = patchi;
387  }
388  else if (procPatchi != -1)
389  {
390  FatalErrorInFunction
391  << "Processor patches should be at end of patch list."
392  << endl
393  << "Have processor patch " << procPatchi
394  << " followed by non-processor patch " << patchi
395  << " in patches " << patches.names()
396  << abort(FatalError);
397  }
398  }
399 
400  return nonEmptyPatchi;
401 }
402 
403 
404 Foam::tmp<Foam::surfaceScalarField> Foam::fvMeshDistribute::generateTestField
405 (
406  const fvMesh& mesh
407 )
408 {
409  const vector testNormal = normalised(vector::one);
410 
411  tmp<surfaceScalarField> tfld
412  (
413  new surfaceScalarField
414  (
415  IOobject
416  (
417  "myFlux",
418  mesh.time().timeName(),
419  mesh,
420  IOobject::NO_READ,
421  IOobject::NO_WRITE
422  ),
423  mesh,
424  dimensionedScalar(dimless, Zero)
425  )
426  );
427  surfaceScalarField& fld = tfld.ref();
428 
429  const surfaceVectorField n(mesh.Sf()/mesh.magSf());
430 
431  forAll(fld, facei)
432  {
433  fld[facei] = (n[facei] & testNormal);
434  }
435 
436  surfaceScalarField::Boundary& fluxBf = fld.boundaryFieldRef();
437  const surfaceVectorField::Boundary& nBf = n.boundaryField();
438 
439  forAll(fluxBf, patchi)
440  {
441  fvsPatchScalarField& fvp = fluxBf[patchi];
442 
443  scalarField newPfld(fvp.size());
444  forAll(newPfld, i)
445  {
446  newPfld[i] = (nBf[patchi][i] & testNormal);
447  }
448  fvp == newPfld;
449  }
450 
451  return tfld;
452 }
453 
454 
455 void Foam::fvMeshDistribute::testField(const surfaceScalarField& fld)
456 {
457  const fvMesh& mesh = fld.mesh();
458 
459  const vector testNormal = normalised(vector::one);
460 
461  const surfaceVectorField n(mesh.Sf()/mesh.magSf());
462 
463  forAll(fld, facei)
464  {
465  scalar cos = (n[facei] & testNormal);
466 
467  if (mag(cos - fld[facei]) > 1e-6)
468  {
469  //FatalErrorInFunction
470  WarningInFunction
471  << "On internal face " << facei << " at "
472  << mesh.faceCentres()[facei]
473  << " the field value is " << fld[facei]
474  << " whereas cos angle of " << testNormal
475  << " with mesh normal " << n[facei]
476  << " is " << cos
477  //<< exit(FatalError);
478  << endl;
479  }
480  }
481  forAll(fld.boundaryField(), patchi)
482  {
483  const fvsPatchScalarField& fvp = fld.boundaryField()[patchi];
484  const fvsPatchVectorField& np = n.boundaryField()[patchi];
485 
486  forAll(fvp, i)
487  {
488  scalar cos = (np[i] & testNormal);
489 
490  if (mag(cos - fvp[i]) > 1e-6)
491  {
492  label facei = fvp.patch().start()+i;
493  //FatalErrorInFunction
494  WarningInFunction
495  << "On face " << facei
496  << " on patch " << fvp.patch().name()
497  << " at " << mesh.faceCentres()[facei]
498  << " the field value is " << fvp[i]
499  << " whereas cos angle of " << testNormal
500  << " with mesh normal " << np[i]
501  << " is " << cos
502  //<< exit(FatalError);
503  << endl;
504  }
505  }
506  }
507 }
508 
509 
510 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::deleteProcPatches
511 (
512  const label destinationPatch
513 )
514 {
515  // Delete all processor patches. Move any processor faces into the last
516  // non-processor patch.
517 
518  // New patchID per boundary faces to be repatched. Is -1 (no change)
519  // or new patchID
520  labelList newPatchID(mesh_.nBoundaryFaces(), -1);
521 
522  for (const polyPatch& pp : mesh_.boundaryMesh())
523  {
524  if (isA<processorPolyPatch>(pp))
525  {
526  if (debug)
527  {
528  Pout<< "Moving all faces of patch " << pp.name()
529  << " into patch " << destinationPatch
530  << endl;
531  }
532 
533  SubList<label>
534  (
535  newPatchID,
536  pp.size(),
537  pp.offset()
538  ) = destinationPatch;
539  }
540  }
541 
542  // Note: order of boundary faces been kept the same since the
543  // destinationPatch is at the end and we have visited the patches in
544  // incremental order.
545  labelListList dummyFaceMaps;
546  autoPtr<mapPolyMesh> map = repatch(newPatchID, dummyFaceMaps);
547 
548 
549  // Delete (now empty) processor patches.
550  {
551  labelList oldToNew(identity(mesh_.boundaryMesh().size()));
552  label newi = 0;
553  // Non processor patches first
554  forAll(mesh_.boundaryMesh(), patchi)
555  {
556  if (!isA<processorPolyPatch>(mesh_.boundaryMesh()[patchi]))
557  {
558  oldToNew[patchi] = newi++;
559  }
560  }
561  label nNonProcPatches = newi;
562 
563  // Processor patches as last
564  forAll(mesh_.boundaryMesh(), patchi)
565  {
566  if (isA<processorPolyPatch>(mesh_.boundaryMesh()[patchi]))
567  {
568  oldToNew[patchi] = newi++;
569  }
570  }
571  fvMeshTools::reorderPatches(mesh_, oldToNew, nNonProcPatches, false);
572  }
573 
574  return map;
575 }
576 
577 
578 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::repatch
579 (
580  const labelList& newPatchID, // per boundary face -1 or new patchID
581  labelListList& constructFaceMap
582 )
583 {
584  polyTopoChange meshMod(mesh_);
585 
586  forAll(newPatchID, bFacei)
587  {
588  if (newPatchID[bFacei] != -1)
589  {
590  label facei = mesh_.nInternalFaces() + bFacei;
591 
592  label zoneID = mesh_.faceZones().whichZone(facei);
593  bool zoneFlip = false;
594 
595  if (zoneID >= 0)
596  {
597  const faceZone& fZone = mesh_.faceZones()[zoneID];
598  zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
599  }
600 
601  meshMod.setAction
602  (
603  polyModifyFace
604  (
605  mesh_.faces()[facei], // modified face
606  facei, // label of face
607  mesh_.faceOwner()[facei], // owner
608  -1, // neighbour
609  false, // face flip
610  newPatchID[bFacei], // patch for face
611  false, // remove from zone
612  zoneID, // zone for face
613  zoneFlip // face flip in zone
614  )
615  );
616  }
617  }
618 
619 
620  // Do mapping of fields from one patchField to the other ourselves since
621  // is currently not supported by updateMesh.
622 
623  // Store boundary fields (we only do this for surfaceFields)
624  PtrList<FieldField<fvsPatchField, scalar>> sFlds;
625  saveBoundaryFields<scalar, surfaceMesh>(sFlds);
626  PtrList<FieldField<fvsPatchField, vector>> vFlds;
627  saveBoundaryFields<vector, surfaceMesh>(vFlds);
628  PtrList<FieldField<fvsPatchField, sphericalTensor>> sptFlds;
629  saveBoundaryFields<sphericalTensor, surfaceMesh>(sptFlds);
630  PtrList<FieldField<fvsPatchField, symmTensor>> sytFlds;
631  saveBoundaryFields<symmTensor, surfaceMesh>(sytFlds);
632  PtrList<FieldField<fvsPatchField, tensor>> tFlds;
633  saveBoundaryFields<tensor, surfaceMesh>(tFlds);
634 
635  // Change the mesh (no inflation). Note: parallel comms allowed.
636  //
637  // NOTE: there is one very particular problem with this ordering.
638  // We first create the processor patches and use these to merge out
639  // shared points (see mergeSharedPoints below). So temporarily points
640  // and edges do not match!
641 
642  // TBD: temporarily unset mesh moving to avoid problems in meshflux
643  // mapping. To be investigated.
644  const bool oldMoving = mesh_.moving(false);
645  autoPtr<mapPolyMesh> mapPtr = meshMod.changeMesh(mesh_, false, true);
646  mesh_.moving(oldMoving);
647  mapPolyMesh& map = *mapPtr;
648 
649 
650  // Update fields. No inflation, parallel sync.
651  mesh_.updateMesh(map);
652 
653  // Map patch fields using stored boundary fields. Note: assumes order
654  // of fields has not changed in object registry!
655  mapBoundaryFields<scalar, surfaceMesh>(map, sFlds);
656  mapBoundaryFields<vector, surfaceMesh>(map, vFlds);
657  mapBoundaryFields<sphericalTensor, surfaceMesh>(map, sptFlds);
658  mapBoundaryFields<symmTensor, surfaceMesh>(map, sytFlds);
659  mapBoundaryFields<tensor, surfaceMesh>(map, tFlds);
660 
661 
662  // Move mesh (since morphing does not do this)
663  if (map.hasMotionPoints())
664  {
665  mesh_.movePoints(map.preMotionPoints());
666  }
667 
668  // Adapt constructMaps.
669 
670  if (debug)
671  {
672  label index = map.reverseFaceMap().find(-1);
673 
674  if (index != -1)
675  {
676  FatalErrorInFunction
677  << "reverseFaceMap contains -1 at index:"
678  << index << endl
679  << "This means that the repatch operation was not just"
680  << " a shuffle?" << abort(FatalError);
681  }
682  }
683 
684  forAll(constructFaceMap, proci)
685  {
686  inplaceRenumberWithFlip
687  (
688  map.reverseFaceMap(),
689  false,
690  true,
691  constructFaceMap[proci]
692  );
693  }
694 
695  return mapPtr;
696 }
697 
698 
699 // Detect shared points. Need processor patches to be present.
700 // Background: when adding bits of mesh one can get points which
701 // share the same position but are only detectable to be topologically
702 // the same point when doing parallel analysis. This routine will
703 // merge those points.
704 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
705 (
706  const labelList& pointToGlobalMaster,
707  labelListList& constructPointMap
708 )
709 {
710  // Find out which sets of points get merged and create a map from
711  // mesh point to unique point.
712 
713  label nShared = 0;
714  forAll(pointToGlobalMaster, pointi)
715  {
716  if (pointToGlobalMaster[pointi] != -1)
717  {
718  nShared++;
719  }
720  }
721 
722  if (debug)
723  {
724  Pout<< "mergeSharedPoints : found " << nShared
725  << " points on processor boundaries" << nl << endl;
726  }
727 
728  Map<label> globalMasterToLocalMaster(2*nShared);
729  Map<label> pointToMaster(2*nShared);
730  label nMatch = 0;
731 
732  forAll(pointToGlobalMaster, pointi)
733  {
734  label globali = pointToGlobalMaster[pointi];
735  if (globali != -1)
736  {
737  const auto iter = globalMasterToLocalMaster.cfind(globali);
738 
739  if (iter.good())
740  {
741  // Matched to existing master
742  pointToMaster.insert(pointi, *iter);
743  nMatch++;
744  }
745  else
746  {
747  // Found first point. Designate as master
748  globalMasterToLocalMaster.insert(globali, pointi);
749  pointToMaster.insert(pointi, pointi);
750  }
751  }
752  }
753 
754  reduce(nMatch, sumOp<label>());
755 
756  if (debug)
757  {
758  Pout<< "mergeSharedPoints : found "
759  << nMatch << " mergeable points" << nl << endl;
760  }
761 
762 
763  if (nMatch == 0)
764  {
765  return nullptr;
766  }
767 
768 
769  polyTopoChange meshMod(mesh_);
770 
771  fvMeshAdder::mergePoints(mesh_, pointToMaster, meshMod);
772 
773  // Change the mesh (no inflation). Note: parallel comms allowed.
774 
775  // TBD: temporarily unset mesh moving to avoid problems in meshflux
776  // mapping. To be investigated.
777  const bool oldMoving = mesh_.moving(false);
778  autoPtr<mapPolyMesh> mapPtr = meshMod.changeMesh(mesh_, false, true);
779  mesh_.moving(oldMoving);
780  mapPolyMesh& map = *mapPtr;
781 
782  // Update fields. No inflation, parallel sync.
783  mesh_.updateMesh(map);
784 
785  // Adapt constructMaps for merged points.
786  forAll(constructPointMap, proci)
787  {
788  labelList& constructMap = constructPointMap[proci];
789 
790  forAll(constructMap, i)
791  {
792  label oldPointi = constructMap[i];
793 
794  label newPointi = map.reversePointMap()[oldPointi];
795 
796  if (newPointi < -1)
797  {
798  constructMap[i] = -newPointi-2;
799  }
800  else if (newPointi >= 0)
801  {
802  constructMap[i] = newPointi;
803  }
804  else
805  {
806  FatalErrorInFunction
807  << "Problem. oldPointi:" << oldPointi
808  << " newPointi:" << newPointi << abort(FatalError);
809  }
810  }
811  }
812 
813  return mapPtr;
814 }
815 
816 
817 void Foam::fvMeshDistribute::getCouplingData
818 (
819  const labelList& distribution,
820  labelList& sourceFace,
821  labelList& sourceProc,
822  labelList& sourcePatch,
823  labelList& sourceNewNbrProc,
824  labelList& sourcePointMaster
825 ) const
826 {
827  // Construct the coupling information for all (boundary) faces and
828  // points
829 
830  const label nBnd = mesh_.nBoundaryFaces();
831  sourceFace.setSize(nBnd);
832  sourceProc.setSize(nBnd);
833  sourcePatch.setSize(nBnd);
834  sourceNewNbrProc.setSize(nBnd);
835 
836  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
837 
838  // Get neighbouring meshFace labels and new processor of coupled boundaries.
839  labelList nbrFaces(nBnd, -1);
840  labelList nbrNewNbrProc(nBnd, -1);
841 
842  forAll(patches, patchi)
843  {
844  const polyPatch& pp = patches[patchi];
845 
846  if (pp.coupled())
847  {
848  label offset = pp.start() - mesh_.nInternalFaces();
849 
850  // Mesh labels of faces on this side
851  forAll(pp, i)
852  {
853  label bndI = offset + i;
854  nbrFaces[bndI] = pp.start()+i;
855  }
856 
857  // Which processor they will end up on
858  SubList<label>(nbrNewNbrProc, pp.size(), offset) =
859  labelUIndList(distribution, pp.faceCells())();
860  }
861  }
862 
863 
864  // Exchange the boundary data
865  syncTools::swapBoundaryFaceList(mesh_, nbrFaces);
866  syncTools::swapBoundaryFaceList(mesh_, nbrNewNbrProc);
867 
868 
869  forAll(patches, patchi)
870  {
871  const polyPatch& pp = patches[patchi];
872  label offset = pp.start() - mesh_.nInternalFaces();
873 
874  if (isA<processorPolyPatch>(pp))
875  {
876  const processorPolyPatch& procPatch =
877  refCast<const processorPolyPatch>(pp);
878 
879  // Check which of the two faces we store.
880 
881  if (procPatch.owner())
882  {
883  // Use my local face labels
884  forAll(pp, i)
885  {
886  label bndI = offset + i;
887  sourceFace[bndI] = pp.start()+i;
888  sourceProc[bndI] = Pstream::myProcNo();
889  sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
890  }
891  }
892  else
893  {
894  // Use my neighbours face labels
895  forAll(pp, i)
896  {
897  label bndI = offset + i;
898  sourceFace[bndI] = nbrFaces[bndI];
899  sourceProc[bndI] = procPatch.neighbProcNo();
900  sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
901  }
902  }
903 
904 
905  label patchi = -1;
906  if (isA<processorCyclicPolyPatch>(pp))
907  {
908  patchi = refCast<const processorCyclicPolyPatch>
909  (
910  pp
911  ).referPatchID();
912  }
913 
914  forAll(pp, i)
915  {
916  label bndI = offset + i;
917  sourcePatch[bndI] = patchi;
918  }
919  }
920  else if (isA<cyclicPolyPatch>(pp))
921  {
922  const cyclicPolyPatch& cpp = refCast<const cyclicPolyPatch>(pp);
923 
924  if (cpp.owner())
925  {
926  forAll(pp, i)
927  {
928  label bndI = offset + i;
929  sourceFace[bndI] = pp.start()+i;
930  sourceProc[bndI] = Pstream::myProcNo();
931  sourcePatch[bndI] = patchi;
932  sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
933  }
934  }
935  else
936  {
937  forAll(pp, i)
938  {
939  label bndI = offset + i;
940  sourceFace[bndI] = nbrFaces[bndI];
941  sourceProc[bndI] = Pstream::myProcNo();
942  sourcePatch[bndI] = patchi;
943  sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
944  }
945  }
946  }
947  else
948  {
949  // Normal (physical) boundary
950  forAll(pp, i)
951  {
952  label bndI = offset + i;
953  sourceFace[bndI] = -1;
954  sourceProc[bndI] = -1;
955  sourcePatch[bndI] = patchi;
956  sourceNewNbrProc[bndI] = -1;
957  }
958  }
959  }
960 
961 
962  // Collect coupled (collocated) points
963  sourcePointMaster.setSize(mesh_.nPoints());
964  sourcePointMaster = -1;
965  {
966  // Assign global master point
967  const globalIndex globalPoints(mesh_.nPoints());
968 
969  const globalMeshData& gmd = mesh_.globalData();
970  const indirectPrimitivePatch& cpp = gmd.coupledPatch();
971  const labelList& meshPoints = cpp.meshPoints();
972  const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
973  const labelListList& slaves = gmd.globalCoPointSlaves();
974 
975  labelList elems(slavesMap.constructSize(), -1);
976  forAll(meshPoints, pointi)
977  {
978  const labelList& slots = slaves[pointi];
979 
980  if (slots.size())
981  {
982  // pointi is a master. Assign a unique label.
983 
984  label globalPointi = globalPoints.toGlobal(meshPoints[pointi]);
985  elems[pointi] = globalPointi;
986  forAll(slots, i)
987  {
988  label sloti = slots[i];
989  if (sloti >= meshPoints.size())
990  {
991  // Filter out local collocated points. We don't want
992  // to merge these
993  elems[slots[i]] = globalPointi;
994  }
995  }
996  }
997  }
998 
999  // Push slave-slot data back to slaves
1000  slavesMap.reverseDistribute(elems.size(), elems, false);
1001 
1002  // Extract back onto mesh
1003  forAll(meshPoints, pointi)
1004  {
1005  sourcePointMaster[meshPoints[pointi]] = elems[pointi];
1006  }
1007  }
1008 }
1009 
1010 
1011 // Subset the neighbourCell/neighbourProc fields
1012 void Foam::fvMeshDistribute::subsetCouplingData
1013 (
1014  const fvMesh& mesh,
1015  const labelList& pointMap,
1016  const labelList& faceMap,
1017  const labelList& cellMap,
1018 
1019  const labelList& oldDistribution,
1020  const labelList& oldFaceOwner,
1021  const labelList& oldFaceNeighbour,
1022  const label oldInternalFaces,
1023 
1024  const labelList& sourceFace,
1025  const labelList& sourceProc,
1026  const labelList& sourcePatch,
1027  const labelList& sourceNewNbrProc,
1028  const labelList& sourcePointMaster,
1029 
1030  labelList& subFace,
1031  labelList& subProc,
1032  labelList& subPatch,
1033  labelList& subNewNbrProc,
1034  labelList& subPointMaster
1035 )
1036 {
1037  subFace.setSize(mesh.nBoundaryFaces());
1038  subProc.setSize(mesh.nBoundaryFaces());
1039  subPatch.setSize(mesh.nBoundaryFaces());
1040  subNewNbrProc.setSize(mesh.nBoundaryFaces());
1041 
1042  forAll(subFace, newBFacei)
1043  {
1044  label newFacei = newBFacei + mesh.nInternalFaces();
1045 
1046  label oldFacei = faceMap[newFacei];
1047 
1048  // Was oldFacei internal face? If so which side did we get.
1049  if (oldFacei < oldInternalFaces)
1050  {
1051  subFace[newBFacei] = oldFacei;
1052  subProc[newBFacei] = Pstream::myProcNo();
1053  subPatch[newBFacei] = -1;
1054 
1055  label oldOwn = oldFaceOwner[oldFacei];
1056  label oldNei = oldFaceNeighbour[oldFacei];
1057 
1058  if (oldOwn == cellMap[mesh.faceOwner()[newFacei]])
1059  {
1060  // We kept the owner side. Where does the neighbour move to?
1061  subNewNbrProc[newBFacei] = oldDistribution[oldNei];
1062  }
1063  else
1064  {
1065  // We kept the neighbour side.
1066  subNewNbrProc[newBFacei] = oldDistribution[oldOwn];
1067  }
1068  }
1069  else
1070  {
1071  // Was boundary face. Take over boundary information
1072  label oldBFacei = oldFacei - oldInternalFaces;
1073 
1074  subFace[newBFacei] = sourceFace[oldBFacei];
1075  subProc[newBFacei] = sourceProc[oldBFacei];
1076  subPatch[newBFacei] = sourcePatch[oldBFacei];
1077  subNewNbrProc[newBFacei] = sourceNewNbrProc[oldBFacei];
1078  }
1079  }
1080 
1081 
1082  subPointMaster = UIndirectList<label>(sourcePointMaster, pointMap);
1083 }
1084 
1085 
1086 // Find cells on mesh whose faceID/procID match the neighbour cell/proc of
1087 // domainMesh. Store the matching face.
1088 void Foam::fvMeshDistribute::findCouples
1089 (
1090  const primitiveMesh& mesh,
1091  const labelList& sourceFace,
1092  const labelList& sourceProc,
1093  const labelList& sourcePatch,
1094 
1095  const label domain,
1096  const primitiveMesh& domainMesh,
1097  const labelList& domainFace,
1098  const labelList& domainProc,
1099  const labelList& domainPatch,
1100 
1101  labelList& masterCoupledFaces,
1102  labelList& slaveCoupledFaces
1103 )
1104 {
1105  // Store domain neighbour as map so we can easily look for pair
1106  // with same face+proc.
1107  labelPairLookup map(domainFace.size());
1108 
1109  forAll(domainProc, bFacei)
1110  {
1111  if (domainProc[bFacei] != -1 && domainPatch[bFacei] == -1)
1112  {
1113  map.insert
1114  (
1115  labelPair(domainFace[bFacei], domainProc[bFacei]),
1116  bFacei
1117  );
1118  }
1119  }
1120 
1121 
1122  // Try to match mesh data.
1123 
1124  masterCoupledFaces.setSize(domainFace.size());
1125  slaveCoupledFaces.setSize(domainFace.size());
1126  label coupledI = 0;
1127 
1128  forAll(sourceFace, bFacei)
1129  {
1130  if (sourceProc[bFacei] != -1 && sourcePatch[bFacei] == -1)
1131  {
1132  labelPair myData(sourceFace[bFacei], sourceProc[bFacei]);
1133 
1134  const auto iter = map.cfind(myData);
1135 
1136  if (iter.good())
1137  {
1138  label nbrBFacei = *iter;
1139 
1140  masterCoupledFaces[coupledI] = mesh.nInternalFaces() + bFacei;
1141  slaveCoupledFaces[coupledI] =
1142  domainMesh.nInternalFaces()
1143  + nbrBFacei;
1144 
1145  coupledI++;
1146  }
1147  }
1148  }
1149 
1150  masterCoupledFaces.setSize(coupledI);
1151  slaveCoupledFaces.setSize(coupledI);
1152 
1153  if (debug)
1154  {
1155  Pout<< "findCouples : found " << coupledI
1156  << " faces that will be stitched" << nl << endl;
1157  }
1158 }
1159 
1160 
1161 void Foam::fvMeshDistribute::findCouples
1162 (
1163  const UPtrList<polyMesh>& meshes,
1164  const PtrList<labelList>& domainSourceFaces,
1165  const PtrList<labelList>& domainSourceProcs,
1166  const PtrList<labelList>& domainSourcePatchs,
1167 
1168  labelListList& localBoundaryFace,
1169  labelListList& remoteFaceProc,
1170  labelListList& remoteBoundaryFace
1171 )
1172 {
1173  // Collect all matching processor face pairs. These are all the
1174  // faces for which we have both sides
1175 
1176  // Pass 0: count number of inter-processor faces
1177  labelList nProcFaces(meshes.size(), 0);
1178  forAll(meshes, meshi)
1179  {
1180  if (meshes.set(meshi))
1181  {
1182  const labelList& domainProc = domainSourceProcs[meshi];
1183  const labelList& domainPatch = domainSourcePatchs[meshi];
1184 
1185  forAll(domainProc, bFacei)
1186  {
1187  if (domainProc[bFacei] != -1 && domainPatch[bFacei] == -1)
1188  {
1189  nProcFaces[meshi]++;
1190  }
1191  }
1192  }
1193  }
1194 
1195  if (debug)
1196  {
1197  Pout<< "fvMeshDistribute::findCouples : nProcFaces:"
1198  << flatOutput(nProcFaces) << endl;
1199  }
1200 
1201 
1202  // Size
1203  List<DynamicList<label>> dynLocalFace(Pstream::nProcs());
1204  List<DynamicList<label>> dynRemoteProc(Pstream::nProcs());
1205  List<DynamicList<label>> dynRemoteFace(Pstream::nProcs());
1206 
1207  forAll(meshes, meshi)
1208  {
1209  if (meshes.set(meshi))
1210  {
1211  dynLocalFace[meshi].setCapacity(nProcFaces[meshi]);
1212  dynRemoteProc[meshi].setCapacity(nProcFaces[meshi]);
1213  dynRemoteFace[meshi].setCapacity(nProcFaces[meshi]);
1214  }
1215  }
1216 
1217 
1218  // Insert all into big map. Find matches
1219  LabelPairMap<labelPair> map(2*sum(nProcFaces));
1220 
1221  nProcFaces = 0;
1222 
1223  forAll(meshes, meshi)
1224  {
1225  if (meshes.set(meshi))
1226  {
1227  const labelList& domainFace = domainSourceFaces[meshi];
1228  const labelList& domainProc = domainSourceProcs[meshi];
1229  const labelList& domainPatch = domainSourcePatchs[meshi];
1230 
1231  forAll(domainProc, bFacei)
1232  {
1233  if (domainProc[bFacei] != -1 && domainPatch[bFacei] == -1)
1234  {
1235  const labelPair key
1236  (
1237  domainProc[bFacei],
1238  domainFace[bFacei]
1239  );
1240  auto fnd = map.find(key);
1241 
1242  if (!fnd.good())
1243  {
1244  // Insert
1245  map.emplace(key, meshi, bFacei);
1246  }
1247  else
1248  {
1249  // Second occurence. Found match.
1250  const label matchProci = fnd().first();
1251  const label matchFacei = fnd().second();
1252 
1253  dynLocalFace[meshi].append(bFacei);
1254  dynRemoteProc[meshi].append(matchProci);
1255  dynRemoteFace[meshi].append(matchFacei);
1256  nProcFaces[meshi]++;
1257 
1258  dynLocalFace[matchProci].append(matchFacei);
1259  dynRemoteProc[matchProci].append(meshi);
1260  dynRemoteFace[matchProci].append(bFacei);
1261  nProcFaces[matchProci]++;
1262  }
1263  }
1264  }
1265  }
1266  }
1267 
1268  if (debug)
1269  {
1270  Pout<< "fvMeshDistribute::findCouples : stored procFaces:"
1271  << map.size() << endl;
1272  }
1273 
1274  localBoundaryFace.setSize(Pstream::nProcs());
1275  remoteFaceProc.setSize(Pstream::nProcs());
1276  remoteBoundaryFace.setSize(Pstream::nProcs());
1277  forAll(meshes, meshi)
1278  {
1279  if (meshes.set(meshi))
1280  {
1281  localBoundaryFace[meshi] = std::move(dynLocalFace[meshi]);
1282  remoteFaceProc[meshi] = std::move(dynRemoteProc[meshi]);
1283  remoteBoundaryFace[meshi] = std::move(dynRemoteFace[meshi]);
1284  }
1285  }
1286 
1287 
1288  if (debug)
1289  {
1290  Pout<< "fvMeshDistribute::findCouples : found matches:"
1291  << flatOutput(nProcFaces) << endl;
1292  }
1293 }
1294 
1295 
1296 // Map data on boundary faces to new mesh (resulting from adding two meshes)
1297 Foam::labelList Foam::fvMeshDistribute::mapBoundaryData
1298 (
1299  const primitiveMesh& mesh, // mesh after adding
1300  const mapAddedPolyMesh& map,
1301  const labelList& boundaryData0, // on mesh before adding
1302  const label nInternalFaces1,
1303  const labelList& boundaryData1 // on added mesh
1304 )
1305 {
1306  labelList newBoundaryData(mesh.nBoundaryFaces());
1307 
1308  forAll(boundaryData0, oldBFacei)
1309  {
1310  label newFacei = map.oldFaceMap()[oldBFacei + map.nOldInternalFaces()];
1311 
1312  // Face still exists (is necessary?) and still boundary face
1313  if (newFacei >= 0 && newFacei >= mesh.nInternalFaces())
1314  {
1315  newBoundaryData[newFacei - mesh.nInternalFaces()] =
1316  boundaryData0[oldBFacei];
1317  }
1318  }
1319 
1320  forAll(boundaryData1, addedBFacei)
1321  {
1322  label newFacei = map.addedFaceMap()[addedBFacei + nInternalFaces1];
1323 
1324  if (newFacei >= 0 && newFacei >= mesh.nInternalFaces())
1325  {
1326  newBoundaryData[newFacei - mesh.nInternalFaces()] =
1327  boundaryData1[addedBFacei];
1328  }
1329  }
1330 
1331  return newBoundaryData;
1332 }
1333 
1334 
1335 Foam::labelList Foam::fvMeshDistribute::mapPointData
1336 (
1337  const primitiveMesh& mesh, // mesh after adding
1338  const mapAddedPolyMesh& map,
1339  const labelList& boundaryData0, // on mesh before adding
1340  const labelList& boundaryData1 // on added mesh
1341 )
1342 {
1343  labelList newBoundaryData(mesh.nPoints());
1344 
1345  forAll(boundaryData0, oldPointi)
1346  {
1347  label newPointi = map.oldPointMap()[oldPointi];
1348 
1349  // Point still exists (is necessary?)
1350  if (newPointi >= 0)
1351  {
1352  newBoundaryData[newPointi] = boundaryData0[oldPointi];
1353  }
1354  }
1355 
1356  forAll(boundaryData1, addedPointi)
1357  {
1358  label newPointi = map.addedPointMap()[addedPointi];
1359 
1360  if (newPointi >= 0)
1361  {
1362  newBoundaryData[newPointi] = boundaryData1[addedPointi];
1363  }
1364  }
1365 
1366  return newBoundaryData;
1367 }
1368 
1369 
1370 // Remove cells. Add all exposed faces to patch oldInternalPatchi
1371 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::doRemoveCells
1372 (
1373  const labelList& cellsToRemove,
1374  const label oldInternalPatchi
1375 )
1376 {
1377  // Mesh change engine
1378  polyTopoChange meshMod(mesh_);
1379 
1380  // Cell removal topo engine. Do NOT synchronize parallel since
1381  // we are doing a local cell removal.
1382  removeCells cellRemover(mesh_, false);
1383 
1384  // Get all exposed faces
1385  labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
1386 
1387  // Insert the topo changes
1388  cellRemover.setRefinement
1389  (
1390  cellsToRemove,
1391  exposedFaces,
1392  labelList(exposedFaces.size(), oldInternalPatchi), // patch for exposed
1393  // faces.
1394  meshMod
1395  );
1396 
1397 
1399  //tmp<surfaceScalarField> sfld(generateTestField(mesh_));
1400 
1401  // Save internal fields (note: not as DimensionedFields since would
1402  // get mapped)
1403  PtrList<Field<scalar>> sFlds;
1404  saveInternalFields(sFlds);
1405  PtrList<Field<vector>> vFlds;
1406  saveInternalFields(vFlds);
1407  PtrList<Field<sphericalTensor>> sptFlds;
1408  saveInternalFields(sptFlds);
1409  PtrList<Field<symmTensor>> sytFlds;
1410  saveInternalFields(sytFlds);
1411  PtrList<Field<tensor>> tFlds;
1412  saveInternalFields(tFlds);
1413 
1414  // Change the mesh. No inflation. Note: no parallel comms allowed.
1415 
1416  // TBD: temporarily unset mesh moving to avoid problems in meshflux
1417  // mapping. To be investigated.
1418  const bool oldMoving = mesh_.moving(false);
1419  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, false);
1420  mesh_.moving(oldMoving);
1421 
1422  // Update fields
1423  mesh_.updateMesh(map());
1424 
1425 
1426  // Any exposed faces in a surfaceField will not be mapped. Map the value
1427  // of these separately (until there is support in all PatchFields for
1428  // mapping from internal faces ...)
1429 
1430  mapExposedFaces(map(), sFlds);
1431  mapExposedFaces(map(), vFlds);
1432  mapExposedFaces(map(), sptFlds);
1433  mapExposedFaces(map(), sytFlds);
1434  mapExposedFaces(map(), tFlds);
1435 
1436 
1438  //testField(sfld);
1439 
1440 
1441  // Move mesh (since morphing does not do this)
1442  if (map().hasMotionPoints())
1443  {
1444  mesh_.movePoints(map().preMotionPoints());
1445  }
1446 
1447 
1448  return map;
1449 }
1450 
1451 
1452 // Delete and add processor patches. Changes mesh and returns per neighbour proc
1453 // the processor patchID.
1454 void Foam::fvMeshDistribute::addProcPatches
1455 (
1456  const labelList& nbrProc, // processor that neighbour is now on
1457  const labelList& referPatchID, // patchID (or -1) I originated from
1458  List<Map<label>>& procPatchID
1459 )
1460 {
1461  // Now use the neighbourFace/Proc to repatch the mesh. These lists
1462  // contain for all current boundary faces the global patchID (for non-proc
1463  // patch) or the processor.
1464 
1465  // Determine a visit order such that the processor patches get added
1466  // in order of increasing neighbour processor (and for same neighbour
1467  // processor (in case of processor cyclics) in order of increasing
1468  // 'refer' patch)
1469  labelList indices;
1470  sortedOrder(nbrProc, indices, lessProcPatches(nbrProc, referPatchID));
1471 
1472  procPatchID.setSize(Pstream::nProcs());
1473 
1474  forAll(indices, i)
1475  {
1476  label bFacei = indices[i];
1477  label proci = nbrProc[bFacei];
1478 
1479  if (proci != -1 && proci != Pstream::myProcNo())
1480  {
1481  if (!procPatchID[proci].found(referPatchID[bFacei]))
1482  {
1483  // No patch for neighbour yet. Is either a normal processor
1484  // patch or a processorCyclic patch.
1485 
1486  if (referPatchID[bFacei] == -1)
1487  {
1488  // Ordinary processor boundary
1489 
1490  processorPolyPatch pp
1491  (
1492  0, // size
1493  mesh_.nFaces(),
1494  mesh_.boundaryMesh().size(),
1495  mesh_.boundaryMesh(),
1496  Pstream::myProcNo(),
1497  proci
1498  );
1499 
1500  procPatchID[proci].insert
1501  (
1502  referPatchID[bFacei],
1503  fvMeshTools::addPatch
1504  (
1505  mesh_,
1506  pp,
1507  dictionary(), // optional per field patchField
1508  processorFvPatchField<scalar>::typeName,
1509  false // not parallel sync
1510  )
1511  );
1512  }
1513  else
1514  {
1515  const coupledPolyPatch& pcPatch
1516  = refCast<const coupledPolyPatch>
1517  (
1518  mesh_.boundaryMesh()[referPatchID[bFacei]]
1519  );
1520  processorCyclicPolyPatch pp
1521  (
1522  0, // size
1523  mesh_.nFaces(),
1524  mesh_.boundaryMesh().size(),
1525  mesh_.boundaryMesh(),
1526  Pstream::myProcNo(),
1527  proci,
1528  pcPatch.name(),
1529  pcPatch.transform()
1530  );
1531 
1532  procPatchID[proci].insert
1533  (
1534  referPatchID[bFacei],
1535  fvMeshTools::addPatch
1536  (
1537  mesh_,
1538  pp,
1539  dictionary(), // optional per field patchField
1540  processorCyclicFvPatchField<scalar>::typeName,
1541  false // not parallel sync
1542  )
1543  );
1544  }
1545  }
1546  }
1547  }
1548 }
1549 
1550 
1551 // Get boundary faces to be repatched. Is -1 or new patchID
1552 Foam::labelList Foam::fvMeshDistribute::getBoundaryPatch
1553 (
1554  const labelList& nbrProc, // new processor per boundary face
1555  const labelList& referPatchID, // patchID (or -1) I originated from
1556  const List<Map<label>>& procPatchID // per proc the new procPatches
1557 )
1558 {
1559  labelList patchIDs(nbrProc);
1560 
1561  forAll(nbrProc, bFacei)
1562  {
1563  if (nbrProc[bFacei] == Pstream::myProcNo())
1564  {
1565  label origPatchi = referPatchID[bFacei];
1566  patchIDs[bFacei] = origPatchi;
1567  }
1568  else if (nbrProc[bFacei] != -1)
1569  {
1570  label origPatchi = referPatchID[bFacei];
1571  patchIDs[bFacei] = procPatchID[nbrProc[bFacei]][origPatchi];
1572  }
1573  else
1574  {
1575  patchIDs[bFacei] = -1;
1576  }
1577  }
1578  return patchIDs;
1579 }
1580 
1581 
1582 // Send mesh and coupling data.
1583 void Foam::fvMeshDistribute::sendMesh
1584 (
1585  const label domain,
1586  const fvMesh& mesh,
1587 
1588  const wordList& pointZoneNames,
1589  const wordList& faceZoneNames,
1590  const wordList& cellZoneNames,
1591 
1592  const labelList& sourceFace,
1593  const labelList& sourceProc,
1594  const labelList& sourcePatch,
1595  const labelList& sourceNewNbrProc,
1596  const labelList& sourcePointMaster,
1597  Ostream& toDomain
1598 )
1599 {
1600  if (debug)
1601  {
1602  Pout<< "Sending to domain " << domain << nl
1603  << " nPoints:" << mesh.nPoints() << nl
1604  << " nFaces:" << mesh.nFaces() << nl
1605  << " nCells:" << mesh.nCells() << nl
1606  << " nPatches:" << mesh.boundaryMesh().size() << nl
1607  << endl;
1608  }
1609 
1610  // Assume sparse, possibly overlapping point zones. Get contents
1611  // in merged-zone indices.
1612  CompactListList<label> zonePoints;
1613  {
1614  const pointZoneMesh& pointZones = mesh.pointZones();
1615 
1616  labelList rowSizes(pointZoneNames.size(), Zero);
1617 
1618  forAll(pointZoneNames, nameI)
1619  {
1620  label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
1621 
1622  if (myZoneID != -1)
1623  {
1624  rowSizes[nameI] = pointZones[myZoneID].size();
1625  }
1626  }
1627  zonePoints.setSize(rowSizes);
1628 
1629  forAll(pointZoneNames, nameI)
1630  {
1631  label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
1632 
1633  if (myZoneID != -1)
1634  {
1635  zonePoints[nameI].deepCopy(pointZones[myZoneID]);
1636  }
1637  }
1638  }
1639 
1640  // Assume sparse, possibly overlapping face zones
1641  CompactListList<label> zoneFaces;
1642  CompactListList<bool> zoneFaceFlip;
1643  {
1644  const faceZoneMesh& faceZones = mesh.faceZones();
1645 
1646  labelList rowSizes(faceZoneNames.size(), Zero);
1647 
1648  forAll(faceZoneNames, nameI)
1649  {
1650  label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
1651 
1652  if (myZoneID != -1)
1653  {
1654  rowSizes[nameI] = faceZones[myZoneID].size();
1655  }
1656  }
1657 
1658  zoneFaces.setSize(rowSizes);
1659  zoneFaceFlip.setSize(rowSizes);
1660 
1661  forAll(faceZoneNames, nameI)
1662  {
1663  label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
1664 
1665  if (myZoneID != -1)
1666  {
1667  zoneFaces[nameI].deepCopy(faceZones[myZoneID]);
1668  zoneFaceFlip[nameI].deepCopy(faceZones[myZoneID].flipMap());
1669  }
1670  }
1671  }
1672 
1673  // Assume sparse, possibly overlapping cell zones
1674  CompactListList<label> zoneCells;
1675  {
1676  const cellZoneMesh& cellZones = mesh.cellZones();
1677 
1678  labelList rowSizes(cellZoneNames.size(), Zero);
1679 
1680  forAll(cellZoneNames, nameI)
1681  {
1682  label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
1683 
1684  if (myZoneID != -1)
1685  {
1686  rowSizes[nameI] = cellZones[myZoneID].size();
1687  }
1688  }
1689 
1690  zoneCells.setSize(rowSizes);
1691 
1692  forAll(cellZoneNames, nameI)
1693  {
1694  label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
1695 
1696  if (myZoneID != -1)
1697  {
1698  zoneCells[nameI].deepCopy(cellZones[myZoneID]);
1699  }
1700  }
1701  }
1703  //labelList cellZoneID;
1704  //if (hasCellZones)
1705  //{
1706  // cellZoneID.setSize(mesh.nCells());
1707  // cellZoneID = -1;
1708  //
1709  // const cellZoneMesh& cellZones = mesh.cellZones();
1710  //
1711  // forAll(cellZones, zoneI)
1712  // {
1713  // labelUIndList(cellZoneID, cellZones[zoneI]) = zoneI;
1714  // }
1715  //}
1716 
1717  // Send
1718  toDomain
1719  << mesh.points()
1720  // Send as (offsets/values)
1721  << CompactListList<label>::pack<face>(mesh.faces())
1722  << mesh.faceOwner()
1723  << mesh.faceNeighbour()
1724  << mesh.boundaryMesh()
1725 
1726  << zonePoints
1727  << zoneFaces
1728  << zoneFaceFlip
1729  << zoneCells
1730 
1731  << sourceFace
1732  << sourceProc
1733  << sourcePatch
1734  << sourceNewNbrProc
1735  << sourcePointMaster;
1736 
1737 
1738  if (debug)
1739  {
1740  Pout<< "Started sending mesh to domain " << domain
1741  << endl;
1742  }
1743 }
1744 
1745 
1746 // Receive mesh. Opposite of sendMesh
1747 Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
1748 (
1749  const label domain,
1750  const wordList& pointZoneNames,
1751  const wordList& faceZoneNames,
1752  const wordList& cellZoneNames,
1753  const Time& runTime,
1754  labelList& domainSourceFace,
1755  labelList& domainSourceProc,
1756  labelList& domainSourcePatch,
1757  labelList& domainSourceNewNbrProc,
1758  labelList& domainSourcePointMaster,
1759  Istream& fromNbr
1760 )
1761 {
1762  pointField domainPoints(fromNbr);
1763  // Receive as (offsets/values), unpack to faceList
1764  faceList domainFaces = CompactListList<label>(fromNbr).unpack<face>();
1765  labelList domainAllOwner(fromNbr);
1766  labelList domainAllNeighbour(fromNbr);
1767  PtrList<entry> patchEntries(fromNbr);
1768 
1769  CompactListList<label> zonePoints(fromNbr);
1770  CompactListList<label> zoneFaces(fromNbr);
1771  CompactListList<bool> zoneFaceFlip(fromNbr);
1772  CompactListList<label> zoneCells(fromNbr);
1773 
1774  fromNbr
1775  >> domainSourceFace
1776  >> domainSourceProc
1777  >> domainSourcePatch
1778  >> domainSourceNewNbrProc
1779  >> domainSourcePointMaster;
1780 
1781  // Construct fvMesh
1782  auto domainMeshPtr = autoPtr<fvMesh>::New
1783  (
1784  IOobject
1785  (
1786  polyMesh::defaultRegion,
1787  runTime.timeName(),
1788  runTime,
1789  IOobject::NO_READ
1790  ),
1791  std::move(domainPoints),
1792  std::move(domainFaces),
1793  std::move(domainAllOwner),
1794  std::move(domainAllNeighbour),
1795  false // no parallel comms
1796  );
1797  fvMesh& domainMesh = *domainMeshPtr;
1798 
1799  polyPatchList patches(patchEntries.size());
1800 
1801  forAll(patchEntries, patchi)
1802  {
1803  patches.set
1804  (
1805  patchi,
1806  polyPatch::New
1807  (
1808  patchEntries[patchi].keyword(),
1809  patchEntries[patchi].dict(),
1810  patchi,
1811  domainMesh.boundaryMesh()
1812  )
1813  );
1814  }
1815  // Add patches; no parallel comms
1816  domainMesh.addFvPatches(patches, false);
1817 
1818  // Construct zones
1819  List<pointZone*> pZonePtrs(pointZoneNames.size());
1820  forAll(pZonePtrs, i)
1821  {
1822  pZonePtrs[i] = new pointZone
1823  (
1824  pointZoneNames[i],
1825  zonePoints[i],
1826  i,
1827  domainMesh.pointZones()
1828  );
1829  }
1830 
1831  List<faceZone*> fZonePtrs(faceZoneNames.size());
1832  forAll(fZonePtrs, i)
1833  {
1834  fZonePtrs[i] = new faceZone
1835  (
1836  faceZoneNames[i],
1837  zoneFaces[i],
1838  zoneFaceFlip[i],
1839  i,
1840  domainMesh.faceZones()
1841  );
1842  }
1843 
1844  List<cellZone*> cZonePtrs(cellZoneNames.size());
1845  forAll(cZonePtrs, i)
1846  {
1847  cZonePtrs[i] = new cellZone
1848  (
1849  cellZoneNames[i],
1850  zoneCells[i],
1851  i,
1852  domainMesh.cellZones()
1853  );
1854  }
1855  domainMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
1856 
1857  return domainMeshPtr;
1858 }
1859 
1860 
1861 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1862 
1863 Foam::fvMeshDistribute::fvMeshDistribute(fvMesh& mesh)
1864 :
1865  mesh_(mesh)
1866 {}
1867 
1868 
1869 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1870 
1871 Foam::labelList Foam::fvMeshDistribute::countCells
1872 (
1873  const labelList& distribution
1874 )
1875 {
1876  labelList nCells(Pstream::nProcs(), Zero);
1877  forAll(distribution, celli)
1878  {
1879  label newProc = distribution[celli];
1880 
1881  if (newProc < 0 || newProc >= Pstream::nProcs())
1882  {
1883  FatalErrorInFunction
1884  << "Distribution should be in range 0.." << Pstream::nProcs()-1
1885  << endl
1886  << "At index " << celli << " distribution:" << newProc
1887  << abort(FatalError);
1888  }
1889  nCells[newProc]++;
1890  }
1891  return nCells;
1892 }
1893 
1894 
1895 Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
1896 (
1897  const labelList& distribution
1898 )
1899 {
1900  // Some checks on distribution
1901  if (distribution.size() != mesh_.nCells())
1902  {
1903  FatalErrorInFunction
1904  << "Size of distribution:"
1905  << distribution.size() << " mesh nCells:" << mesh_.nCells()
1906  << abort(FatalError);
1907  }
1908 
1909 
1910  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1911 
1912  // Check all processors have same non-proc patches in same order.
1913  if (patches.checkParallelSync(true))
1914  {
1915  FatalErrorInFunction
1916  << "This application requires all non-processor patches"
1917  << " to be present in the same order on all patches" << nl
1918  << "followed by the processor patches (which of course are unique)."
1919  << nl
1920  << "Local patches:" << mesh_.boundaryMesh().names()
1921  << abort(FatalError);
1922  }
1923 
1924  // Save some data for mapping later on
1925  const label nOldPoints(mesh_.nPoints());
1926  const label nOldFaces(mesh_.nFaces());
1927  const label nOldCells(mesh_.nCells());
1928  labelList oldPatchStarts(patches.size());
1929  labelList oldPatchNMeshPoints(patches.size());
1930  forAll(patches, patchi)
1931  {
1932  oldPatchStarts[patchi] = patches[patchi].start();
1933  oldPatchNMeshPoints[patchi] = patches[patchi].nPoints();
1934  }
1935 
1936 
1937 
1938  // Short circuit trivial case.
1939  if (!Pstream::parRun())
1940  {
1941  // Collect all maps and return
1942  return autoPtr<mapDistributePolyMesh>::New
1943  (
1944  mesh_,
1945 
1946  nOldPoints,
1947  nOldFaces,
1948  nOldCells,
1949  std::move(oldPatchStarts),
1950  std::move(oldPatchNMeshPoints),
1951 
1952  labelListList(one(), identity(mesh_.nPoints())), //subPointMap
1953  labelListList(one(), identity(mesh_.nFaces())), //subFaceMap
1954  labelListList(one(), identity(mesh_.nCells())), //subCellMap
1955  labelListList(one(), identity(patches.size())), //subPatchMap
1956 
1957  labelListList(one(), identity(mesh_.nPoints())), //pointMap
1958  labelListList(one(), identity(mesh_.nFaces())), //faceMap
1959  labelListList(one(), identity(mesh_.nCells())), //cellMap
1960  labelListList(one(), identity(patches.size())) //patchMap
1961  );
1962  }
1963 
1964 
1965  // Collect any zone names over all processors and shuffle zones accordingly
1966  // Note that this is not necessary for redistributePar since that already
1967  // checks for it. However other use (e.g. mesh generators) might not.
1968  const wordList pointZoneNames(mergeWordList(mesh_.pointZones().names()));
1969  reorderZones<pointZone>(pointZoneNames, mesh_.pointZones());
1970 
1971  const wordList faceZoneNames(mergeWordList(mesh_.faceZones().names()));
1972  reorderZones<faceZone>(faceZoneNames, mesh_.faceZones());
1973 
1974  const wordList cellZoneNames(mergeWordList(mesh_.cellZones().names()));
1975  reorderZones<cellZone>(cellZoneNames, mesh_.cellZones());
1976 
1977 
1978  // Local environment of all boundary faces
1979  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1980 
1981  // A face is uniquely defined by
1982  // - proc
1983  // - local face no
1984  //
1985  // To glue the parts of meshes which can get sent from anywhere we
1986  // need to know on boundary faces what the above tuple on both sides is.
1987  // So we need to maintain:
1988  // - original face
1989  // - original processor id (= trivial)
1990  // For coupled boundaries (where the faces are 'duplicate') we take the
1991  // lowest numbered processor as the data to store.
1992  //
1993  // Additionally to create the procboundaries we need to know where the owner
1994  // cell on the other side now is: newNeighbourProc.
1995  //
1996 
1997  // physical boundary:
1998  // sourceProc = -1
1999  // sourceNewNbrProc = -1
2000  // sourceFace = -1
2001  // sourcePatch = patchID
2002  // processor boundary:
2003  // sourceProc = proc (on owner side)
2004  // sourceNewNbrProc = distribution of coupled cell
2005  // sourceFace = face (on owner side)
2006  // sourcePatch = -1
2007  // ?cyclic:
2008  // ? sourceProc = proc
2009  // ? sourceNewNbrProc = distribution of coupled cell
2010  // ? sourceFace = face (on owner side)
2011  // ? sourcePatch = patchID
2012  // processor-cyclic boundary:
2013  // sourceProc = proc (on owner side)
2014  // sourceNewNbrProc = distribution of coupled cell
2015  // sourceFace = face (on owner side)
2016  // sourcePatch = patchID
2017 
2018  labelList sourcePatch;
2019  labelList sourceFace;
2020  labelList sourceProc;
2021  labelList sourceNewNbrProc;
2022  labelList sourcePointMaster;
2023  getCouplingData
2024  (
2025  distribution,
2026  sourceFace,
2027  sourceProc,
2028  sourcePatch,
2029  sourceNewNbrProc,
2030  sourcePointMaster
2031  );
2032 
2033 
2034  // Remove meshPhi. Since this would otherwise disappear anyway
2035  // during topo changes and we have to guarantee that all the fields
2036  // can be sent.
2037  mesh_.clearOut();
2038  mesh_.resetMotion();
2039 
2040  // Get data to send. Make sure is synchronised
2041 
2042  HashTable<wordList> allFieldNames;
2043 
2044  getFieldNames<volScalarField>(mesh_, allFieldNames);
2045  getFieldNames<volVectorField>(mesh_, allFieldNames);
2046  getFieldNames<volSphericalTensorField>(mesh_, allFieldNames);
2047  getFieldNames<volSymmTensorField>(mesh_, allFieldNames);
2048  getFieldNames<volTensorField>(mesh_, allFieldNames);
2049 
2050  getFieldNames<surfaceScalarField>(mesh_, allFieldNames);
2051  getFieldNames<surfaceVectorField>(mesh_, allFieldNames);
2052  getFieldNames<surfaceSphericalTensorField>(mesh_, allFieldNames);
2053  getFieldNames<surfaceSymmTensorField>(mesh_, allFieldNames);
2054  getFieldNames<surfaceTensorField>(mesh_, allFieldNames);
2055 
2056  getFieldNames<volScalarField::Internal>
2057  (
2058  mesh_,
2059  allFieldNames,
2060  volScalarField::typeName
2061  );
2062  getFieldNames<volVectorField::Internal>
2063  (
2064  mesh_,
2065  allFieldNames,
2066  volVectorField::typeName
2067  );
2068  getFieldNames<volSphericalTensorField::Internal>
2069  (
2070  mesh_,
2071  allFieldNames,
2072  volSphericalTensorField::typeName
2073  );
2074  getFieldNames<volSymmTensorField::Internal>
2075  (
2076  mesh_,
2077  allFieldNames,
2078  volSymmTensorField::typeName
2079  );
2080  getFieldNames<volTensorField::Internal>
2081  (
2082  mesh_,
2083  allFieldNames,
2084  volTensorField::typeName
2085  );
2086 
2087 
2088  // Find patch to temporarily put exposed and processor faces into.
2089  const label oldInternalPatchi = findNonEmptyPatch();
2090 
2091 
2092  // Delete processor patches, starting from the back. Move all faces into
2093  // oldInternalPatchi.
2094  labelList repatchFaceMap;
2095  {
2096  autoPtr<mapPolyMesh> repatchMap = deleteProcPatches(oldInternalPatchi);
2097 
2098  // Store face map (only face ordering that changed)
2099  repatchFaceMap = repatchMap().faceMap();
2100 
2101  // Reorder all boundary face data (sourceProc, sourceFace etc.)
2102  labelList bFaceMap
2103  (
2104  SubList<label>
2105  (
2106  repatchMap().reverseFaceMap(),
2107  mesh_.nBoundaryFaces(),
2108  mesh_.nInternalFaces()
2109  )
2110  - mesh_.nInternalFaces()
2111  );
2112 
2113  inplaceReorder(bFaceMap, sourceFace);
2114  inplaceReorder(bFaceMap, sourceProc);
2115  inplaceReorder(bFaceMap, sourcePatch);
2116  inplaceReorder(bFaceMap, sourceNewNbrProc);
2117  }
2118 
2119 
2120 
2121  // Print a bit.
2122  if (debug)
2123  {
2124  Pout<< nl << "MESH WITH PROC PATCHES DELETED:" << endl;
2125  printMeshInfo(mesh_);
2126  printFieldInfo<volScalarField>(mesh_);
2127  printFieldInfo<volVectorField>(mesh_);
2128  printFieldInfo<volSphericalTensorField>(mesh_);
2129  printFieldInfo<volSymmTensorField>(mesh_);
2130  printFieldInfo<volTensorField>(mesh_);
2131  printFieldInfo<surfaceScalarField>(mesh_);
2132  printFieldInfo<surfaceVectorField>(mesh_);
2133  printFieldInfo<surfaceSphericalTensorField>(mesh_);
2134  printFieldInfo<surfaceSymmTensorField>(mesh_);
2135  printFieldInfo<surfaceTensorField>(mesh_);
2136  printIntFieldInfo<volScalarField::Internal>(mesh_);
2137  printIntFieldInfo<volVectorField::Internal>(mesh_);
2138  printIntFieldInfo<volSphericalTensorField::Internal>(mesh_);
2139  printIntFieldInfo<volSymmTensorField::Internal>(mesh_);
2140  printIntFieldInfo<volTensorField::Internal>(mesh_);
2141  Pout<< nl << endl;
2142  }
2143 
2144 
2145 
2146  // Maps from subsetted mesh (that is sent) back to original maps
2147  labelListList subCellMap(Pstream::nProcs());
2148  labelListList subFaceMap(Pstream::nProcs());
2149  labelListList subPointMap(Pstream::nProcs());
2150  labelListList subPatchMap(Pstream::nProcs());
2151  // Maps from subsetted mesh to reconstructed mesh
2152  labelListList constructCellMap(Pstream::nProcs());
2153  labelListList constructFaceMap(Pstream::nProcs());
2154  labelListList constructPointMap(Pstream::nProcs());
2155  labelListList constructPatchMap(Pstream::nProcs());
2156 
2157 
2158  // Find out schedule
2159  // ~~~~~~~~~~~~~~~~~
2160 
2161  labelList nSendCells(countCells(distribution));
2162  labelList nRecvCells(Pstream::nProcs());
2163  UPstream::allToAll(nSendCells, nRecvCells);
2164 
2165  // Allocate buffers
2166  PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
2167 
2168 
2169  // What to send to neighbouring domains
2170  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2171 
2172  // Disable parallel.
2173  const bool oldParRun = UPstream::parRun(false);
2174 
2175  forAll(nSendCells, recvProc)
2176  {
2177  if (recvProc != Pstream::myProcNo() && nSendCells[recvProc] > 0)
2178  {
2179  // Send to recvProc
2180 
2181  if (debug)
2182  {
2183  Pout<< nl
2184  << "SUBSETTING FOR DOMAIN " << recvProc
2185  << " cells to send:"
2186  << nSendCells[recvProc]
2187  << nl << endl;
2188  }
2189 
2190  // Pstream for sending mesh and fields
2191  //OPstream str(Pstream::commsTypes::blocking, recvProc);
2192  UOPstream str(recvProc, pBufs);
2193 
2194  // Mesh subsetting engine - subset the cells of the current domain.
2195  fvMeshSubset subsetter
2196  (
2197  mesh_,
2198  recvProc,
2199  distribution,
2200  oldInternalPatchi, // oldInternalFaces patch
2201  false // no parallel sync
2202  );
2203 
2204  subCellMap[recvProc] = subsetter.cellMap();
2205  subFaceMap[recvProc] = subsetter.faceFlipMap();
2206  inplaceRenumberWithFlip
2207  (
2208  repatchFaceMap,
2209  false, // oldToNew has flip
2210  true, // subFaceMap has flip
2211  subFaceMap[recvProc]
2212  );
2213  subPointMap[recvProc] = subsetter.pointMap();
2214  subPatchMap[recvProc] = subsetter.patchMap();
2215 
2216 
2217  // Subset the boundary fields (owner/neighbour/processor)
2218  labelList procSourceFace;
2219  labelList procSourceProc;
2220  labelList procSourcePatch;
2221  labelList procSourceNewNbrProc;
2222  labelList procSourcePointMaster;
2223 
2224  subsetCouplingData
2225  (
2226  subsetter.subMesh(),
2227  subsetter.pointMap(), // from subMesh to mesh
2228  subsetter.faceMap(), // ,, ,,
2229  subsetter.cellMap(), // ,, ,,
2230 
2231  distribution, // old mesh distribution
2232  mesh_.faceOwner(), // old owner
2233  mesh_.faceNeighbour(),
2234  mesh_.nInternalFaces(),
2235 
2236  sourceFace,
2237  sourceProc,
2238  sourcePatch,
2239  sourceNewNbrProc,
2240  sourcePointMaster,
2241 
2242  procSourceFace,
2243  procSourceProc,
2244  procSourcePatch,
2245  procSourceNewNbrProc,
2246  procSourcePointMaster
2247  );
2248 
2249 
2250  // Send to neighbour
2251  sendMesh
2252  (
2253  recvProc,
2254  subsetter.subMesh(),
2255 
2256  pointZoneNames,
2257  faceZoneNames,
2258  cellZoneNames,
2259 
2260  procSourceFace,
2261  procSourceProc,
2262  procSourcePatch,
2263  procSourceNewNbrProc,
2264  procSourcePointMaster,
2265 
2266  str
2267  );
2268 
2269  // volFields
2270  sendFields<volScalarField>
2271  (
2272  recvProc,
2273  allFieldNames,
2274  subsetter,
2275  str
2276  );
2277  sendFields<volVectorField>
2278  (
2279  recvProc,
2280  allFieldNames,
2281  subsetter,
2282  str
2283  );
2284  sendFields<volSphericalTensorField>
2285  (
2286  recvProc,
2287  allFieldNames,
2288  subsetter,
2289  str
2290  );
2291  sendFields<volSymmTensorField>
2292  (
2293  recvProc,
2294  allFieldNames,
2295  subsetter,
2296  str
2297  );
2298  sendFields<volTensorField>
2299  (
2300  recvProc,
2301  allFieldNames,
2302  subsetter,
2303  str
2304  );
2305 
2306  // surfaceFields
2307  sendFields<surfaceScalarField>
2308  (
2309  recvProc,
2310  allFieldNames,
2311  subsetter,
2312  str
2313  );
2314  sendFields<surfaceVectorField>
2315  (
2316  recvProc,
2317  allFieldNames,
2318  subsetter,
2319  str
2320  );
2321  sendFields<surfaceSphericalTensorField>
2322  (
2323  recvProc,
2324  allFieldNames,
2325  subsetter,
2326  str
2327  );
2328  sendFields<surfaceSymmTensorField>
2329  (
2330  recvProc,
2331  allFieldNames,
2332  subsetter,
2333  str
2334  );
2335  sendFields<surfaceTensorField>
2336  (
2337  recvProc,
2338  allFieldNames,
2339  subsetter,
2340  str
2341  );
2342 
2343  // Dimensioned fields
2344  sendFields<volScalarField::Internal>
2345  (
2346  recvProc,
2347  allFieldNames,
2348  subsetter,
2349  str
2350  );
2351  sendFields<volVectorField::Internal>
2352  (
2353  recvProc,
2354  allFieldNames,
2355  subsetter,
2356  str
2357  );
2358  sendFields<volSphericalTensorField::Internal>
2359  (
2360  recvProc,
2361  allFieldNames,
2362  subsetter,
2363  str
2364  );
2365  sendFields<volSymmTensorField::Internal>
2366  (
2367  recvProc,
2368  allFieldNames,
2369  subsetter,
2370  str
2371  );
2372  sendFields<volTensorField::Internal>
2373  (
2374  recvProc,
2375  allFieldNames,
2376  subsetter,
2377  str
2378  );
2379  }
2380  }
2381 
2382 
2383  UPstream::parRun(oldParRun); // Restore parallel state
2384 
2385  if (debug)
2386  {
2387  Pout<< "Starting sending" << endl;
2388  }
2389 
2390  pBufs.finishedSends();
2391 
2392  if (debug)
2393  {
2394  Pout<< "Finished sending and receiving : "
2395  << flatOutput(pBufs.recvDataCounts()) << endl;
2396  }
2397 
2398 
2399  // Subset the part that stays
2400  // ~~~~~~~~~~~~~~~~~~~~~~~~~~
2401 
2402  {
2403  // Save old mesh maps before changing mesh
2404  const labelList oldFaceOwner(mesh_.faceOwner());
2405  const labelList oldFaceNeighbour(mesh_.faceNeighbour());
2406  const label oldInternalFaces = mesh_.nInternalFaces();
2407 
2408  // Remove cells.
2409  autoPtr<mapPolyMesh> subMap
2410  (
2411  doRemoveCells
2412  (
2413  select(false, distribution, Pstream::myProcNo()),
2414  oldInternalPatchi
2415  )
2416  );
2417 
2418  // Addressing from subsetted mesh
2419  subCellMap[Pstream::myProcNo()] = subMap().cellMap();
2420  subFaceMap[Pstream::myProcNo()] = renumber
2421  (
2422  repatchFaceMap,
2423  subMap().faceMap()
2424  );
2425  // Insert the sign bit from face flipping
2426  labelList& faceMap = subFaceMap[Pstream::myProcNo()];
2427  forAll(faceMap, faceI)
2428  {
2429  faceMap[faceI] += 1;
2430  }
2431  const labelHashSet& flip = subMap().flipFaceFlux();
2432  for (const label facei : flip)
2433  {
2434  faceMap[facei] = -faceMap[facei];
2435  }
2436  subPointMap[Pstream::myProcNo()] = subMap().pointMap();
2437  subPatchMap[Pstream::myProcNo()] = identity(patches.size());
2438 
2439  // Subset the mesh data: neighbourCell/neighbourProc fields
2440  labelList domainSourceFace;
2441  labelList domainSourceProc;
2442  labelList domainSourcePatch;
2443  labelList domainSourceNewNbrProc;
2444  labelList domainSourcePointMaster;
2445 
2446  subsetCouplingData
2447  (
2448  mesh_, // new mesh
2449  subMap().pointMap(), // from new to original mesh
2450  subMap().faceMap(), // from new to original mesh
2451  subMap().cellMap(),
2452 
2453  distribution, // distribution before subsetting
2454  oldFaceOwner, // owner before subsetting
2455  oldFaceNeighbour, // neighbour ,,
2456  oldInternalFaces, // nInternalFaces ,,
2457 
2458  sourceFace,
2459  sourceProc,
2460  sourcePatch,
2461  sourceNewNbrProc,
2462  sourcePointMaster,
2463 
2464  domainSourceFace,
2465  domainSourceProc,
2466  domainSourcePatch,
2467  domainSourceNewNbrProc,
2468  domainSourcePointMaster
2469  );
2470 
2471  sourceFace.transfer(domainSourceFace);
2472  sourceProc.transfer(domainSourceProc);
2473  sourcePatch.transfer(domainSourcePatch);
2474  sourceNewNbrProc.transfer(domainSourceNewNbrProc);
2475  sourcePointMaster.transfer(domainSourcePointMaster);
2476  }
2477 
2478 
2479  // Print a bit.
2480  if (debug)
2481  {
2482  Pout<< nl << "STARTING MESH:" << endl;
2483  printMeshInfo(mesh_);
2484  printFieldInfo<volScalarField>(mesh_);
2485  printFieldInfo<volVectorField>(mesh_);
2486  printFieldInfo<volSphericalTensorField>(mesh_);
2487  printFieldInfo<volSymmTensorField>(mesh_);
2488  printFieldInfo<volTensorField>(mesh_);
2489  printFieldInfo<surfaceScalarField>(mesh_);
2490  printFieldInfo<surfaceVectorField>(mesh_);
2491  printFieldInfo<surfaceSphericalTensorField>(mesh_);
2492  printFieldInfo<surfaceSymmTensorField>(mesh_);
2493  printFieldInfo<surfaceTensorField>(mesh_);
2494  printIntFieldInfo<volScalarField::Internal>(mesh_);
2495  printIntFieldInfo<volVectorField::Internal>(mesh_);
2496  printIntFieldInfo<volSphericalTensorField::Internal>(mesh_);
2497  printIntFieldInfo<volSymmTensorField::Internal>(mesh_);
2498  printIntFieldInfo<volTensorField::Internal>(mesh_);
2499  Pout<< nl << endl;
2500  }
2501 
2502 
2503 
2504  // Receive and add what was sent
2505  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2506 
2507  // Disable parallel. Original state already known.
2508  UPstream::parRun(false);
2509 
2510  PtrList<labelList> domainSourceFaces(Pstream::nProcs());
2511  PtrList<labelList> domainSourceProcs(Pstream::nProcs());
2512  PtrList<labelList> domainSourcePatchs(Pstream::nProcs());
2513  PtrList<labelList> domainSourceNewNbrProcs(Pstream::nProcs());
2514  PtrList<labelList> domainSourcePointMasters(Pstream::nProcs());
2515 
2516  PtrList<fvMesh> domainMeshPtrs(Pstream::nProcs());
2517 
2518  PtrList<PtrList<volScalarField>> vsfs(Pstream::nProcs());
2519  PtrList<PtrList<volVectorField>> vvfs(Pstream::nProcs());
2520  PtrList<PtrList<volSphericalTensorField>> vsptfs(Pstream::nProcs());
2521  PtrList<PtrList<volSymmTensorField>> vsytfs(Pstream::nProcs());
2522  PtrList<PtrList<volTensorField>> vtfs(Pstream::nProcs());
2523 
2524  PtrList<PtrList<surfaceScalarField>> ssfs(Pstream::nProcs());
2525  PtrList<PtrList<surfaceVectorField>> svfs(Pstream::nProcs());
2526  PtrList<PtrList<surfaceSphericalTensorField>> ssptfs
2527  (
2528  Pstream::nProcs()
2529  );
2530  PtrList<PtrList<surfaceSymmTensorField>> ssytfs(Pstream::nProcs());
2531  PtrList<PtrList<surfaceTensorField>> stfs(Pstream::nProcs());
2532 
2533  PtrList<PtrList<volScalarField::Internal>> dsfs(Pstream::nProcs());
2534  PtrList<PtrList<volVectorField::Internal>> dvfs(Pstream::nProcs());
2535  PtrList<PtrList<volSphericalTensorField::Internal>> dstfs
2536  (
2537  Pstream::nProcs()
2538  );
2539  PtrList<PtrList<volSymmTensorField::Internal>> dsytfs
2540  (
2541  Pstream::nProcs()
2542  );
2543  PtrList<PtrList<volTensorField::Internal>> dtfs(Pstream::nProcs());
2544 
2545  forAll(nRecvCells, sendProc)
2546  {
2547  // Did processor sendProc send anything to me?
2548  if (sendProc != Pstream::myProcNo() && nRecvCells[sendProc] > 0)
2549  {
2550  if (debug)
2551  {
2552  Pout<< nl
2553  << "RECEIVING FROM DOMAIN " << sendProc
2554  << " cells to receive:"
2555  << nRecvCells[sendProc]
2556  << nl << endl;
2557  }
2558 
2559 
2560  // Pstream for receiving mesh and fields
2561  UIPstream str(sendProc, pBufs);
2562 
2563 
2564  // Receive from sendProc
2565  domainSourceFaces.set(sendProc, new labelList(0));
2566  labelList& domainSourceFace = domainSourceFaces[sendProc];
2567 
2568  domainSourceProcs.set(sendProc, new labelList(0));
2569  labelList& domainSourceProc = domainSourceProcs[sendProc];
2570 
2571  domainSourcePatchs.set(sendProc, new labelList(0));
2572  labelList& domainSourcePatch = domainSourcePatchs[sendProc];
2573 
2574  domainSourceNewNbrProcs.set(sendProc, new labelList(0));
2575  labelList& domainSourceNewNbrProc =
2576  domainSourceNewNbrProcs[sendProc];
2577 
2578  domainSourcePointMasters.set(sendProc, new labelList(0));
2579  labelList& domainSourcePointMaster =
2580  domainSourcePointMasters[sendProc];
2581 
2582  // Opposite of sendMesh
2583  {
2584  autoPtr<fvMesh> domainMeshPtr = receiveMesh
2585  (
2586  sendProc,
2587  pointZoneNames,
2588  faceZoneNames,
2589  cellZoneNames,
2590 
2591  const_cast<Time&>(mesh_.time()),
2592  domainSourceFace,
2593  domainSourceProc,
2594  domainSourcePatch,
2595  domainSourceNewNbrProc,
2596  domainSourcePointMaster,
2597  str
2598  );
2599  domainMeshPtrs.set(sendProc, domainMeshPtr.ptr());
2600  fvMesh& domainMesh = domainMeshPtrs[sendProc];
2601  // Force construction of various on mesh.
2602  //(void)domainMesh.globalData();
2603 
2604 
2605  // Receive fields. Read as single dictionary because
2606  // of problems reading consecutive fields from single stream.
2607  dictionary fieldDicts(str);
2608 
2609  // Vol fields
2610  vsfs.set(sendProc, new PtrList<volScalarField>(0));
2611  receiveFields<volScalarField>
2612  (
2613  sendProc,
2614  allFieldNames,
2615  domainMesh,
2616  vsfs[sendProc],
2617  fieldDicts
2618  );
2619  vvfs.set(sendProc, new PtrList<volVectorField>(0));
2620  receiveFields<volVectorField>
2621  (
2622  sendProc,
2623  allFieldNames,
2624  domainMesh,
2625  vvfs[sendProc],
2626  fieldDicts
2627  );
2628  vsptfs.set
2629  (
2630  sendProc,
2631  new PtrList<volSphericalTensorField>(0)
2632  );
2633  receiveFields<volSphericalTensorField>
2634  (
2635  sendProc,
2636  allFieldNames,
2637  domainMesh,
2638  vsptfs[sendProc],
2639  fieldDicts
2640  );
2641  vsytfs.set(sendProc, new PtrList<volSymmTensorField>(0));
2642  receiveFields<volSymmTensorField>
2643  (
2644  sendProc,
2645  allFieldNames,
2646  domainMesh,
2647  vsytfs[sendProc],
2648  fieldDicts
2649  );
2650  vtfs.set(sendProc, new PtrList<volTensorField>(0));
2651  receiveFields<volTensorField>
2652  (
2653  sendProc,
2654  allFieldNames,
2655  domainMesh,
2656  vtfs[sendProc],
2657  fieldDicts
2658  );
2659 
2660  // Surface fields
2661  ssfs.set(sendProc, new PtrList<surfaceScalarField>(0));
2662  receiveFields<surfaceScalarField>
2663  (
2664  sendProc,
2665  allFieldNames,
2666  domainMesh,
2667  ssfs[sendProc],
2668  fieldDicts
2669  );
2670  svfs.set(sendProc, new PtrList<surfaceVectorField>(0));
2671  receiveFields<surfaceVectorField>
2672  (
2673  sendProc,
2674  allFieldNames,
2675  domainMesh,
2676  svfs[sendProc],
2677  fieldDicts
2678  );
2679  ssptfs.set
2680  (
2681  sendProc,
2682  new PtrList<surfaceSphericalTensorField>(0)
2683  );
2684  receiveFields<surfaceSphericalTensorField>
2685  (
2686  sendProc,
2687  allFieldNames,
2688  domainMesh,
2689  ssptfs[sendProc],
2690  fieldDicts
2691  );
2692  ssytfs.set(sendProc, new PtrList<surfaceSymmTensorField>(0));
2693  receiveFields<surfaceSymmTensorField>
2694  (
2695  sendProc,
2696  allFieldNames,
2697  domainMesh,
2698  ssytfs[sendProc],
2699  fieldDicts
2700  );
2701  stfs.set(sendProc, new PtrList<surfaceTensorField>(0));
2702  receiveFields<surfaceTensorField>
2703  (
2704  sendProc,
2705  allFieldNames,
2706  domainMesh,
2707  stfs[sendProc],
2708  fieldDicts
2709  );
2710 
2711  // Dimensioned fields
2712  dsfs.set
2713  (
2714  sendProc,
2715  new PtrList<volScalarField::Internal>(0)
2716  );
2717  receiveFields<volScalarField::Internal>
2718  (
2719  sendProc,
2720  allFieldNames,
2721  domainMesh,
2722  dsfs[sendProc],
2723  fieldDicts
2724  );
2725  dvfs.set
2726  (
2727  sendProc,
2728  new PtrList<volVectorField::Internal>(0)
2729  );
2730  receiveFields<volVectorField::Internal>
2731  (
2732  sendProc,
2733  allFieldNames,
2734  domainMesh,
2735  dvfs[sendProc],
2736  fieldDicts
2737  );
2738  dstfs.set
2739  (
2740  sendProc,
2741  new PtrList<volSphericalTensorField::Internal>(0)
2742  );
2743  receiveFields<volSphericalTensorField::Internal>
2744  (
2745  sendProc,
2746  allFieldNames,
2747  domainMesh,
2748  dstfs[sendProc],
2749  fieldDicts
2750  );
2751  dsytfs.set
2752  (
2753  sendProc,
2754  new PtrList<volSymmTensorField::Internal>(0)
2755  );
2756  receiveFields<volSymmTensorField::Internal>
2757  (
2758  sendProc,
2759  allFieldNames,
2760  domainMesh,
2761  dsytfs[sendProc],
2762  fieldDicts
2763  );
2764  dtfs.set
2765  (
2766  sendProc,
2767  new PtrList<volTensorField::Internal>(0)
2768  );
2769  receiveFields<volTensorField::Internal>
2770  (
2771  sendProc,
2772  allFieldNames,
2773  domainMesh,
2774  dtfs[sendProc],
2775  fieldDicts
2776  );
2777  }
2778  }
2779  }
2780 
2781  // Clear out storage
2782  pBufs.clear();
2783 
2784 
2785  // Set up pointers to meshes so we can include our mesh_
2786  UPtrList<polyMesh> meshes(domainMeshPtrs.size());
2787  UPtrList<fvMesh> fvMeshes(domainMeshPtrs.size());
2788  forAll(domainMeshPtrs, proci)
2789  {
2790  if (domainMeshPtrs.set(proci))
2791  {
2792  meshes.set(proci, &domainMeshPtrs[proci]);
2793  fvMeshes.set(proci, &domainMeshPtrs[proci]);
2794  }
2795  }
2796 
2797  // 'Receive' from myself
2798  {
2799  meshes.set(Pstream::myProcNo(), &mesh_);
2800  fvMeshes.set(Pstream::myProcNo(), &mesh_);
2801 
2802  //domainSourceFaces.set(Pstream::myProcNo(), std::move(sourceFace));
2803  domainSourceFaces.set(Pstream::myProcNo(), new labelList(0));
2804  domainSourceFaces[Pstream::myProcNo()] = sourceFace;
2805 
2806  domainSourceProcs.set(Pstream::myProcNo(), new labelList(0));
2807  //std::move(sourceProc));
2808  domainSourceProcs[Pstream::myProcNo()] = sourceProc;
2809 
2810  domainSourcePatchs.set(Pstream::myProcNo(), new labelList(0));
2811  //, std::move(sourcePatch));
2812  domainSourcePatchs[Pstream::myProcNo()] = sourcePatch;
2813 
2814  domainSourceNewNbrProcs.set(Pstream::myProcNo(), new labelList(0));
2815  domainSourceNewNbrProcs[Pstream::myProcNo()] = sourceNewNbrProc;
2816 
2817  domainSourcePointMasters.set(Pstream::myProcNo(), new labelList(0));
2818  domainSourcePointMasters[Pstream::myProcNo()] = sourcePointMaster;
2819  }
2820 
2821 
2822  // Find matching faces that need to be stitched
2823  labelListList localBoundaryFace(Pstream::nProcs());
2824  labelListList remoteFaceProc(Pstream::nProcs());
2825  labelListList remoteBoundaryFace(Pstream::nProcs());
2826  findCouples
2827  (
2828  meshes,
2829  domainSourceFaces,
2830  domainSourceProcs,
2831  domainSourcePatchs,
2832 
2833  localBoundaryFace,
2834  remoteFaceProc,
2835  remoteBoundaryFace
2836  );
2837 
2838 
2839  const label nOldInternalFaces = mesh_.nInternalFaces();
2840  const labelList oldFaceOwner(mesh_.faceOwner());
2841 
2842  // TBD: temporarily unset mesh moving to avoid problems in meshflux
2843  // mapping. To be investigated.
2844  const bool oldMoving = mesh_.moving(false);
2845 
2846  fvMeshAdder::add
2847  (
2848  Pstream::myProcNo(), // index of mesh to modify (== mesh_)
2849  fvMeshes,
2850  oldFaceOwner,
2851 
2852  // Coupling info
2853  localBoundaryFace,
2854  remoteFaceProc,
2855  remoteBoundaryFace,
2856 
2857  constructPatchMap,
2858  constructCellMap,
2859  constructFaceMap,
2860  constructPointMap
2861  );
2862 
2863  mesh_.moving(oldMoving);
2864 
2865 
2866  if (debug)
2867  {
2868  Pout<< nl << "ADDED REMOTE MESHES:" << endl;
2869  printMeshInfo(mesh_);
2870  printFieldInfo<volScalarField>(mesh_);
2871  printFieldInfo<volVectorField>(mesh_);
2872  printFieldInfo<volSphericalTensorField>(mesh_);
2873  printFieldInfo<volSymmTensorField>(mesh_);
2874  printFieldInfo<volTensorField>(mesh_);
2875  printFieldInfo<surfaceScalarField>(mesh_);
2876  printFieldInfo<surfaceVectorField>(mesh_);
2877  printFieldInfo<surfaceSphericalTensorField>(mesh_);
2878  printFieldInfo<surfaceSymmTensorField>(mesh_);
2879  printFieldInfo<surfaceTensorField>(mesh_);
2880  printIntFieldInfo<volScalarField::Internal>(mesh_);
2881  printIntFieldInfo<volVectorField::Internal>(mesh_);
2882  printIntFieldInfo<volSphericalTensorField::Internal>(mesh_);
2883  printIntFieldInfo<volSymmTensorField::Internal>(mesh_);
2884  printIntFieldInfo<volTensorField::Internal>(mesh_);
2885  Pout<< nl << endl;
2886  }
2887 
2888  {
2889  //- Combine sourceProc, sourcePatch, sourceFace
2890  sourceProc.setSize(mesh_.nBoundaryFaces());
2891  sourceProc = -1;
2892  sourcePatch.setSize(mesh_.nBoundaryFaces());
2893  sourcePatch = -1;
2894  sourceFace.setSize(mesh_.nBoundaryFaces());
2895  sourceFace = -1;
2896  sourceNewNbrProc.setSize(mesh_.nBoundaryFaces());
2897  sourceNewNbrProc = -1;
2898  sourcePointMaster.setSize(mesh_.nPoints());
2899  sourcePointMaster = -1;
2900 
2901  if (mesh_.nPoints() > 0)
2902  {
2903  forAll(meshes, meshi)
2904  {
2905  if (domainSourceFaces.set(meshi))
2906  {
2907  const label nIntFaces =
2908  (
2909  meshi == Pstream::myProcNo()
2910  ? nOldInternalFaces
2911  : meshes[meshi].nInternalFaces()
2912  );
2913  const labelList& faceOwner
2914  (
2915  meshi == Pstream::myProcNo()
2916  ? oldFaceOwner
2917  : meshes[meshi].faceOwner()
2918  );
2919 
2920  labelList& faceMap = constructFaceMap[meshi];
2921  const labelList& cellMap = constructCellMap[meshi];
2922 
2923  const labelList& domainSourceFace =
2924  domainSourceFaces[meshi];
2925  const labelList& domainSourceProc =
2926  domainSourceProcs[meshi];
2927  const labelList& domainSourcePatch =
2928  domainSourcePatchs[meshi];
2929  const labelList& domainSourceNewNbr =
2930  domainSourceNewNbrProcs[meshi];
2931  UIndirectList<label>
2932  (
2933  sourcePointMaster,
2934  constructPointMap[meshi]
2935  ) = domainSourcePointMasters[meshi];
2936 
2937 
2938  forAll(domainSourceFace, bFacei)
2939  {
2940  const label oldFacei = bFacei+nIntFaces;
2941  const label allFacei = faceMap[oldFacei];
2942  const label allbFacei = allFacei-mesh_.nInternalFaces();
2943 
2944  if (allbFacei >= 0)
2945  {
2946  sourceProc[allbFacei] = domainSourceProc[bFacei];
2947  sourcePatch[allbFacei] = domainSourcePatch[bFacei];
2948  sourceFace[allbFacei] = domainSourceFace[bFacei];
2949  sourceNewNbrProc[allbFacei] =
2950  domainSourceNewNbr[bFacei];
2951  }
2952  }
2953 
2954 
2955  // Add flip to constructFaceMap
2956  forAll(faceMap, oldFacei)
2957  {
2958  const label allFacei = faceMap[oldFacei];
2959  const label allOwn = mesh_.faceOwner()[allFacei];
2960 
2961  if (cellMap[faceOwner[oldFacei]] == allOwn)
2962  {
2963  // Master face
2964  faceMap[oldFacei] += 1;
2965  }
2966  else
2967  {
2968  // Slave face. Flip.
2969  faceMap[oldFacei] = -faceMap[oldFacei] - 1;
2970  }
2971  }
2972  }
2973  }
2974  }
2975  }
2976 
2977 
2978  UPstream::parRun(oldParRun); // Restore parallel state
2979 
2980 
2981  // Print a bit.
2982  if (debug)
2983  {
2984  Pout<< nl << "REDISTRIBUTED MESH:" << endl;
2985  printMeshInfo(mesh_);
2986  printFieldInfo<volScalarField>(mesh_);
2987  printFieldInfo<volVectorField>(mesh_);
2988  printFieldInfo<volSphericalTensorField>(mesh_);
2989  printFieldInfo<volSymmTensorField>(mesh_);
2990  printFieldInfo<volTensorField>(mesh_);
2991  printFieldInfo<surfaceScalarField>(mesh_);
2992  printFieldInfo<surfaceVectorField>(mesh_);
2993  printFieldInfo<surfaceSphericalTensorField>(mesh_);
2994  printFieldInfo<surfaceSymmTensorField>(mesh_);
2995  printFieldInfo<surfaceTensorField>(mesh_);
2996  printIntFieldInfo<volScalarField::Internal>(mesh_);
2997  printIntFieldInfo<volVectorField::Internal>(mesh_);
2998  printIntFieldInfo<volSphericalTensorField::Internal>(mesh_);
2999  printIntFieldInfo<volSymmTensorField::Internal>(mesh_);
3000  printIntFieldInfo<volTensorField::Internal>(mesh_);
3001  Pout<< nl << endl;
3002  }
3003 
3004 
3005  // See if any originally shared points need to be merged. Note: does
3006  // parallel comms. After this points and edges should again be consistent.
3007  mergeSharedPoints(sourcePointMaster, constructPointMap);
3008 
3009 
3010  // Add processorPatches
3011  // ~~~~~~~~~~~~~~~~~~~~
3012 
3013  // Per neighbour processor, per originating patch, the patchID
3014  // For faces resulting from internal faces or normal processor patches
3015  // the originating patch is -1. For cyclics this is the cyclic patchID.
3016  List<Map<label>> procPatchID;
3017 
3018  // Add processor and processorCyclic patches.
3019  addProcPatches(sourceNewNbrProc, sourcePatch, procPatchID);
3020 
3021  // Put faces into correct patch. Note that we now have proper
3022  // processorPolyPatches again so repatching will take care of coupled face
3023  // ordering.
3024 
3025  // Get boundary faces to be repatched. Is -1 or new patchID
3026  labelList newPatchID
3027  (
3028  getBoundaryPatch
3029  (
3030  sourceNewNbrProc,
3031  sourcePatch,
3032  procPatchID
3033  )
3034  );
3035 
3036  // Change patches. Since this might change ordering of coupled faces
3037  // we also need to adapt our constructMaps.
3038  repatch(newPatchID, constructFaceMap);
3039 
3040  // Bit of hack: processorFvPatchField does not get reset since created
3041  // from nothing so explicitly reset.
3042  initPatchFields<volScalarField, processorFvPatchField<scalar>>
3043  (
3044  Zero
3045  );
3046  initPatchFields<volVectorField, processorFvPatchField<vector>>
3047  (
3048  Zero
3049  );
3050  initPatchFields
3051  <
3052  volSphericalTensorField,
3053  processorFvPatchField<sphericalTensor>
3054  >
3055  (
3056  Zero
3057  );
3058  initPatchFields<volSymmTensorField, processorFvPatchField<symmTensor>>
3059  (
3060  Zero
3061  );
3062  initPatchFields<volTensorField, processorFvPatchField<tensor>>
3063  (
3064  Zero
3065  );
3066 
3067 
3068  mesh_.setInstance(mesh_.time().timeName());
3069 
3070 
3071  // Print a bit
3072  if (debug)
3073  {
3074  Pout<< nl << "FINAL MESH:" << endl;
3075  printMeshInfo(mesh_);
3076  printFieldInfo<volScalarField>(mesh_);
3077  printFieldInfo<volVectorField>(mesh_);
3078  printFieldInfo<volSphericalTensorField>(mesh_);
3079  printFieldInfo<volSymmTensorField>(mesh_);
3080  printFieldInfo<volTensorField>(mesh_);
3081  printFieldInfo<surfaceScalarField>(mesh_);
3082  printFieldInfo<surfaceVectorField>(mesh_);
3083  printFieldInfo<surfaceSphericalTensorField>(mesh_);
3084  printFieldInfo<surfaceSymmTensorField>(mesh_);
3085  printFieldInfo<surfaceTensorField>(mesh_);
3086  printIntFieldInfo<volScalarField::Internal>(mesh_);
3087  printIntFieldInfo<volVectorField::Internal>(mesh_);
3088  printIntFieldInfo<volSphericalTensorField::Internal>(mesh_);
3089  printIntFieldInfo<volSymmTensorField::Internal>(mesh_);
3090  printIntFieldInfo<volTensorField::Internal>(mesh_);
3091  Pout<< nl << endl;
3092  }
3093 
3094  // Collect all maps and return
3095  return autoPtr<mapDistributePolyMesh>::New
3096  (
3097  mesh_,
3098 
3099  nOldPoints,
3100  nOldFaces,
3101  nOldCells,
3102  std::move(oldPatchStarts),
3103  std::move(oldPatchNMeshPoints),
3104 
3105  std::move(subPointMap),
3106  std::move(subFaceMap),
3107  std::move(subCellMap),
3108  std::move(subPatchMap),
3109 
3110  std::move(constructPointMap),
3111  std::move(constructFaceMap),
3112  std::move(constructCellMap),
3113  std::move(constructPatchMap),
3114 
3115  true, // subFaceMap has flip
3116  true // constructFaceMap has flip
3117  );
3118 }
3119 
3120 
3121 // ************************************************************************* //
surfaceFields.H
Foam::surfaceFields.
Foam::sign
dimensionedScalar sign(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:159
Foam::fvsPatchVectorField
fvsPatchField< vector > fvsPatchVectorField
Definition: fvsPatchFieldsFwd.H:42
Foam::OSstream::name
virtual const fileName & name() const
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
Definition: OSstream.H:128
Foam::fvsPatchScalarField
fvsPatchField< scalar > fvsPatchScalarField
Definition: fvsPatchFieldsFwd.H:39
dict
dictionary dict
Definition: searchingEngine.H:11
Foam::fvMeshTools::reorderPatches
static void reorderPatches(fvMesh &, const labelList &oldToNew, const label nPatches, const bool validBoundary)
Reorder and remove trailing patches.
Definition: fvMeshTools.C:317
Foam::faceZoneMesh
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
Definition: faceZoneMeshFwd.H:38
patchIDs
const labelList patchIDs(pbm.patchSet(polyPatchNames, false, true).sortedToc())
Foam::fvMesh::Sf
const surfaceVectorField & Sf() const
Return cell face area vectors.
Definition: fvMeshGeometry.C:300
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
Definition: DimensionedFieldFunctions.C:295
Foam::PtrListOps::names
List< word > names(const UPtrList< T > &list, const UnaryMatchPredicate &matcher)
List of names generated by calling name() for each list item and filtered for matches.
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:125
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::primitiveMesh::nPoints
label nPoints() const noexcept
Number of mesh points.
Definition: primitiveMeshI.H:30
Foam::FatalError
error FatalError
Error stream (stdout output on all processes), with additional &#39;FOAM FATAL ERROR&#39; header text and sta...
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:578
Foam::polyMesh::faceNeighbour
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1117
Foam::sortedOrder
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
Foam::IOobject::name
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:162
Foam::tmp::ref
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
Definition: tmpI.H:244
Foam::List::set
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:472
Foam::renumber
IntListType renumber(const labelUList &oldToNew, const IntListType &input)
Renumber the values (not the indices) of a list.
Definition: ListOpsTemplates.C:26
Foam::BitOps::select
List< bool > select(const label n, const labelUList &locations)
Construct a selection list of bools (all false) with the given pre-size, subsequently add specified l...
Definition: BitOps.C:134
Foam::nl
constexpr char nl
The newline &#39;\n&#39; character (0x0a)
Definition: Ostream.H:49
Foam::inplaceReorder
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
Definition: ListOpsTemplates.C:113
Foam::polyBoundaryMesh::start
label start() const noexcept
The start label of boundary faces in the polyMesh face list.
Definition: polyBoundaryMesh.C:633
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: IndirectList.H:65
runTime
engineTime & runTime
Definition: createEngineTime.H:13
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:487
Foam::UPstream::parRun
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:1004
Foam::lessProcPatches
Less function class that can be used for sorting processor patches.
Definition: fvMeshDistribute.C:54
Foam::boundBox
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:63
Foam::CompactListList::pack
static CompactListList< T > pack(const UList< SubListType > &lists, const bool checkOverflow=false)
Construct by packing together the list of lists.
Foam::GeometricField< scalar, fvsPatchField, surfaceMesh >::Boundary
GeometricBoundaryField< scalar, fvsPatchField, surfaceMesh > Boundary
Type of boundary fields.
Definition: GeometricField.H:99
Foam::IOobjectOption::NO_WRITE
Ignore writing from objectRegistry::writeObject()
Definition: IOobjectOption.H:86
Foam::normalised
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:673
Foam::dimless
const dimensionSet dimless
Dimensionless.
Definition: dimensionSets.C:182
Foam::UPstream::myProcNo
static int myProcNo(const label communicator=worldComm)
Rank of this process in the communicator (starting from masterNo()). Can be negative if the process i...
Definition: UPstream.H:1029
Foam::fvMesh::time
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:362
Foam::primitiveMesh::nFaces
label nFaces() const noexcept
Number of mesh faces.
Definition: primitiveMeshI.H:83
Foam::indirectPrimitivePatch
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field...
Definition: indirectPrimitivePatch.H:51
Foam::labelListList
List< labelList > labelListList
List of labelList.
Definition: labelList.H:38
Foam::volSphericalTensorField
GeometricField< sphericalTensor, fvPatchField, volMesh > volSphericalTensorField
Definition: volFieldsFwd.H:84
Foam::labelUList
UList< label > labelUList
A UList of labels.
Definition: UList.H:78
Foam::faceMap
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
Definition: blockMeshMergeTopological.C:87
Foam::Pstream::broadcast
static void broadcast(Type &value, const label comm=UPstream::worldComm)
Broadcast content (contiguous or non-contiguous) to all processes in communicator.
Definition: PstreamBroadcast.C:28
ListOps.H
Various functions to operate on Lists.
Foam::polyMesh::points
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1073
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:414
Foam::labelHashSet
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
Foam::FieldBase::typeName
static const char *const typeName
Typename for Field.
Definition: Field.H:86
Foam::faceList
List< face > faceList
List of faces.
Definition: faceListFwd.H:39
Foam::HashTableOps::values
List< T > values(const HashTable< T, Key, Hash > &tbl, const bool doSort=false)
List of values from HashTable, optionally sorted.
Definition: HashOps.H:164
Foam::polyPatchList
PtrList< polyPatch > polyPatchList
Store lists of polyPatch as a PtrList.
Definition: polyPatch.H:56
Foam::HashTable::size
label size() const noexcept
The number of elements in table.
Definition: HashTable.H:331
Foam::UPstream::nProcs
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run. ...
Definition: UPstream.H:1020
Foam::fvMeshDistribute
Sends/receives parts of mesh+fvfields to neighbouring processors. Used in load balancing.
Definition: fvMeshDistribute.H:66
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:38
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. Uses the specified communication schedule...
Definition: PstreamGatherList.C:41
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:258
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:52
Foam::lessProcPatches::operator()
bool operator()(const label a, const label b)
Definition: fvMeshDistribute.C:67
Foam::identity
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:31
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Foam::polyMeshAdder::mergePoints
static void mergePoints(const polyMesh &, const Map< label > &pointToMaster, polyTopoChange &meshMod)
Helper: Merge points.
Definition: polyMeshAdder.C:2191
Foam::polyMesh::boundaryMesh
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
Definition: polyMesh.H:584
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:46
Foam::fvMeshTools::addPatch
static label addPatch(fvMesh &mesh, const polyPatch &patch, const dictionary &patchFieldDict, const word &defaultPatchFieldType, const bool validBoundary)
Add patch. Inserts patch before all processor patches.
Definition: fvMeshTools.C:38
Foam::polyMesh::defaultRegion
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:397
Foam::polyBoundaryMesh::names
wordList names() const
Return a list of patch names.
Definition: polyBoundaryMesh.C:571
Foam::UPtrList::size
label size() const noexcept
The number of entries in the list.
Definition: UPtrListI.H:113
meshes
Foam::PtrList< Foam::fvMesh > meshes(regionNames.size())
Foam::fvMeshAdder::add
static autoPtr< mapAddedPolyMesh > add(fvMesh &mesh0, const fvMesh &mesh1, const faceCoupleInfo &coupleInfo, const bool validBoundary=true, const bool fullyMapped=false)
Inplace add mesh to fvMesh. Maps all stored fields. Returns map.
Definition: fvMeshAdder.C:67
Foam::polyMesh::faceOwner
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1111
Foam::primitiveMesh::nInternalFaces
label nInternalFaces() const noexcept
Number of internal faces.
Definition: primitiveMeshI.H:71
Foam::pointZoneMesh
ZoneMesh< pointZone, polyMesh > pointZoneMesh
A ZoneMesh with the type pointZone.
Definition: pointZoneMeshFwd.H:38
Foam::polyMesh::faces
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1098
Foam::distribution
Accumulating histogram of values. Specified bin resolution automatic generation of bins...
Definition: distribution.H:57
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:139
Foam::PtrList::set
const T * set(const label i) const
Return const pointer to element (can be nullptr), or nullptr for out-of-range access (ie...
Definition: PtrList.H:159
labelPairHashes.H
A HashTable to objects of type <T> with a labelPair key. The hashing is based on labelPair (FixedList...
Foam::fvMesh::magSf
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
Definition: fvMeshGeometry.C:311
Foam::polyBoundaryMesh::checkParallelSync
bool checkParallelSync(const bool report=false) const
Check whether all procs have all patches and in same order.
Definition: polyBoundaryMesh.C:1013
Foam::UPstream::allToAll
static void allToAll(const UList< int32_t > &sendData, UList< int32_t > &recvData, const label communicator=worldComm)
Exchange integer data with all processors (in the communicator).
Definition: UPstreamAllToAll.C:42
Foam::vector
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
Foam::Time::timeName
static word timeName(const scalar t, const int precision=precision_)
Return time name of given scalar time formatted with the given precision.
Definition: Time.C:770
Foam::ensightOutput::debug
int debug
Static debugging option.
Foam::labelPair
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:51
Foam::polyMesh::faceZones
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
Definition: polyMesh.H:646
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
Foam::fvMeshDistribute::distribute
autoPtr< mapDistributePolyMesh > distribute(const labelList &dist)
Send cells to neighbours according to distribution.
Definition: fvMeshDistribute.C:1891
Foam::DimensionedField::typeName
const word typeName("volScalarField::Internal")
fld
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< ' ';}gmvFile<< nl;for(const word &name :lagrangianScalarNames){ IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Definition: gmvOutputLagrangian.H:23
Foam::primitiveMesh::faceCentres
const vectorField & faceCentres() const
Definition: primitiveMeshFaceCentresAndAreas.C:71
Foam::polyMesh::pointZones
const pointZoneMesh & pointZones() const noexcept
Return point zone mesh.
Definition: polyMesh.H:638
Foam::wordList
List< word > wordList
List of word.
Definition: fileName.H:58
Foam::fvMeshDistribute::printMeshInfo
static void printMeshInfo(const fvMesh &)
Print some info on mesh.
Definition: fvMeshDistribute.C:220
Foam::inplaceRenumber
void inplaceRenumber(const labelUList &oldToNew, IntListType &input)
Inplace renumber the values (not the indices) of a list.
Definition: ListOpsTemplates.C:50
Foam::fvMeshDistribute::generateTestField
static tmp< surfaceScalarField > generateTestField(const fvMesh &)
Generate a test field on faces.
Definition: fvMeshDistribute.C:400
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:414
newPointi
label newPointi
Definition: readKivaGrid.H:496
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:33
Foam::primitiveMesh::nCells
label nCells() const noexcept
Number of mesh cells.
Definition: primitiveMeshI.H:89
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:103
Foam::fvBoundaryMesh
Foam::fvBoundaryMesh.
Definition: fvBoundaryMesh.H:53
Foam::fvMeshDistribute::printCoupleInfo
static void printCoupleInfo(const primitiveMesh &, const labelList &, const labelList &, const labelList &, const labelList &)
Print some info on coupling data.
Definition: fvMeshDistribute.C:280
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:79
Foam::UPstream::master
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
Definition: UPstream.H:1037
patches
const polyBoundaryMesh & patches
Definition: convertProcessorPatches.H:54
Foam::IOobjectOption::NO_READ
Nothing to be read.
Definition: IOobjectOption.H:60
Foam::lessProcPatches::lessProcPatches
lessProcPatches(const labelList &nbrProc, const labelList &referPatchID)
Definition: fvMeshDistribute.C:61
Foam::labelPairLookup
HashTable< label, labelPair, Foam::Hash< labelPair > > labelPairLookup
This is a Map of a labelPair to a label. Used for e.g. for face1, face2 to shared edge...
Definition: labelPairHashes.H:63
Foam::polyMesh::cellZones
const cellZoneMesh & cellZones() const noexcept
Return cell zone mesh.
Definition: polyMesh.H:654
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
Definition: PstreamReduceOps.H:49
Foam::UPstream::commsTypes::nonBlocking
"nonBlocking" : (MPI_Isend, MPI_Irecv)
Foam::fvMesh::boundary
const fvBoundaryMesh & boundary() const noexcept
Return reference to boundary mesh.
Definition: fvMesh.H:397
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::polyPatch::New
static autoPtr< polyPatch > New(const word &patchType, const word &name, const label size, const label start, const label index, const polyBoundaryMesh &bm)
Return pointer to a new patch created on freestore from components.
Definition: polyPatchNew.C:28
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
Definition: HashPtrTable.H:48
forAllReverse
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: stdFoam.H:430
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:62
Foam::tmp
A class for managing temporary objects.
Definition: HashPtrTable.H:50
Foam::autoPtr::New
static autoPtr< T > New(Args &&... args)
Construct autoPtr with forwarding arguments.
Definition: autoPtr.H:178
Foam::polyPatch
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:69
Foam::primitiveMesh::nBoundaryFaces
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
Definition: primitiveMeshI.H:77
Foam::fvMeshDistribute::countCells
static labelList countCells(const labelList &)
Helper function: count cells per processor in wanted distribution.
Definition: fvMeshDistribute.C:1867
Foam::fvMeshDistribute::testField
static void testField(const surfaceScalarField &)
Check whether field consistent with face orientation.
Definition: fvMeshDistribute.C:450
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported...
Definition: IOobject.H:171
found
bool found
Definition: TABSMDCalcMethod2.H:32
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::cellZoneMesh
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.
Definition: cellZoneMeshFwd.H:38
Foam::HashTable::set
bool set(const Key &key, const T &obj)
Copy assign a new entry, overwriting existing entries.
Definition: HashTableI.H:174
pp
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:26
Foam::syncTools::swapBoundaryFaceList
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &faceValues)
Swap coupled boundary face values. Uses eqOp.
Definition: syncTools.H:485
Foam::flatOutput
FlatOutput::OutputAdaptor< Container, Delimiters > flatOutput(const Container &obj, Delimiters delim)
Global flatOutput() function with specified output delimiters.
Definition: FlatOutput.H:225
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:133

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.14