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 | Copyright (C) 2011-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "fvMeshDistribute.H"
28 #include "fvMeshAdder.H"
29 #include "faceCoupleInfo.H"
30 #include "processorFvPatchField.H"
31 #include "processorFvsPatchField.H"
34 #include "polyTopoChange.H"
35 #include "removeCells.H"
36 #include "polyModifyFace.H"
37 #include "polyRemovePoint.H"
38 #include "mapDistributePolyMesh.H"
39 #include "surfaceFields.H"
40 #include "syncTools.H"
41 #include "CompactListList.H"
42 #include "fvMeshTools.H"
43 #include "ListOps.H"
44 
45 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
46 
47 namespace Foam
48 {
49  defineTypeNameAndDebug(fvMeshDistribute, 0);
50 
51 //- Less function class that can be used for sorting processor patches
53 {
54  const labelList& nbrProc_;
55  const labelList& referPatchID_;
56 
57 public:
58 
59  lessProcPatches( const labelList& nbrProc, const labelList& referPatchID)
60  :
61  nbrProc_(nbrProc),
62  referPatchID_(referPatchID)
63  {}
64 
65  bool operator()(const label a, const label b)
66  {
67  if (nbrProc_[a] < nbrProc_[b])
68  {
69  return true;
70  }
71  else if (nbrProc_[a] > nbrProc_[b])
72  {
73  return false;
74  }
75  else
76  {
77  // Equal neighbour processor
78  return referPatchID_[a] < referPatchID_[b];
79  }
80  }
81 };
82 
83 }
84 
85 
86 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
87 
88 void Foam::fvMeshDistribute::inplaceRenumberWithFlip
89 (
90  const labelUList& oldToNew,
91  const bool oldToNewHasFlip,
92  const bool lstHasFlip,
93  labelUList& lst
94 )
95 {
96  if (!lstHasFlip && !oldToNewHasFlip)
97  {
98  Foam::inplaceRenumber(oldToNew, lst);
99  }
100  else
101  {
102  // Either input data or map encodes sign so result encodes sign
103 
104  forAll(lst, elemI)
105  {
106  // Extract old value and sign
107  label val = lst[elemI];
108  label sign = 1;
109  if (lstHasFlip)
110  {
111  if (val > 0)
112  {
113  val = val-1;
114  }
115  else if (val < 0)
116  {
117  val = -val-1;
118  sign = -1;
119  }
120  else
121  {
123  << "Problem : zero value " << val
124  << " at index " << elemI << " out of " << lst.size()
125  << " list with flip bit" << exit(FatalError);
126  }
127  }
128 
129 
130  // Lookup new value and possibly change sign
131  label newVal = oldToNew[val];
132 
133  if (oldToNewHasFlip)
134  {
135  if (newVal > 0)
136  {
137  newVal = newVal-1;
138  }
139  else if (newVal < 0)
140  {
141  newVal = -newVal-1;
142  sign = -sign;
143  }
144  else
145  {
147  << "Problem : zero value " << newVal
148  << " at index " << elemI << " out of "
149  << oldToNew.size()
150  << " list with flip bit" << exit(FatalError);
151  }
152  }
153 
154 
155  // Encode new value and sign
156  lst[elemI] = sign*(newVal+1);
157  }
158  }
159 }
160 
161 
162 Foam::labelList Foam::fvMeshDistribute::select
163 (
164  const bool selectEqual,
165  const labelList& values,
166  const label value
167 )
168 {
169  label n = 0;
170 
171  forAll(values, i)
172  {
173  if (selectEqual == (values[i] == value))
174  {
175  n++;
176  }
177  }
178 
179  labelList indices(n);
180  n = 0;
181 
182  forAll(values, i)
183  {
184  if (selectEqual == (values[i] == value))
185  {
186  indices[n++] = i;
187  }
188  }
189  return indices;
190 }
191 
192 
193 // Check all procs have same names and in exactly same order.
194 void Foam::fvMeshDistribute::checkEqualWordList
195 (
196  const string& msg,
197  const wordList& lst
198 )
199 {
200  List<wordList> allNames(Pstream::nProcs());
201  allNames[Pstream::myProcNo()] = lst;
202  Pstream::gatherList(allNames);
203  Pstream::scatterList(allNames);
204 
205  for (label proci = 1; proci < Pstream::nProcs(); proci++)
206  {
207  if (allNames[proci] != allNames[0])
208  {
210  << "When checking for equal " << msg.c_str() << " :" << endl
211  << "processor0 has:" << allNames[0] << endl
212  << "processor" << proci << " has:" << allNames[proci] << endl
213  << msg.c_str() << " need to be synchronised on all processors."
214  << exit(FatalError);
215  }
216  }
217 }
218 
219 
220 Foam::wordList Foam::fvMeshDistribute::mergeWordList(const wordList& procNames)
221 {
222  List<wordList> allNames(Pstream::nProcs());
223  allNames[Pstream::myProcNo()] = procNames;
224  Pstream::gatherList(allNames);
225  Pstream::scatterList(allNames);
226 
227  HashSet<word> mergedNames;
228  forAll(allNames, proci)
229  {
230  forAll(allNames[proci], i)
231  {
232  mergedNames.insert(allNames[proci][i]);
233  }
234  }
235  return mergedNames.toc();
236 }
237 
238 
239 // Print some info on mesh.
241 {
242  Pout<< "Primitives:" << nl
243  << " points :" << mesh.nPoints() << nl
244  << " bb :" << boundBox(mesh.points(), false) << nl
245  << " internalFaces:" << mesh.nInternalFaces() << nl
246  << " faces :" << mesh.nFaces() << nl
247  << " cells :" << mesh.nCells() << nl;
248 
249  const fvBoundaryMesh& patches = mesh.boundary();
250 
251  Pout<< "Patches:" << endl;
252  forAll(patches, patchi)
253  {
254  const polyPatch& pp = patches[patchi].patch();
255 
256  Pout<< " " << patchi << " name:" << pp.name()
257  << " size:" << pp.size()
258  << " start:" << pp.start()
259  << " type:" << pp.type()
260  << endl;
261  }
262 
263  if (mesh.pointZones().size())
264  {
265  Pout<< "PointZones:" << endl;
266  forAll(mesh.pointZones(), zoneI)
267  {
268  const pointZone& pz = mesh.pointZones()[zoneI];
269  Pout<< " " << zoneI << " name:" << pz.name()
270  << " size:" << pz.size()
271  << endl;
272  }
273  }
274  if (mesh.faceZones().size())
275  {
276  Pout<< "FaceZones:" << endl;
277  forAll(mesh.faceZones(), zoneI)
278  {
279  const faceZone& fz = mesh.faceZones()[zoneI];
280  Pout<< " " << zoneI << " name:" << fz.name()
281  << " size:" << fz.size()
282  << endl;
283  }
284  }
285  if (mesh.cellZones().size())
286  {
287  Pout<< "CellZones:" << endl;
288  forAll(mesh.cellZones(), zoneI)
289  {
290  const cellZone& cz = mesh.cellZones()[zoneI];
291  Pout<< " " << zoneI << " name:" << cz.name()
292  << " size:" << cz.size()
293  << endl;
294  }
295  }
296 }
297 
298 
300 (
301  const primitiveMesh& mesh,
302  const labelList& sourceFace,
303  const labelList& sourceProc,
304  const labelList& sourcePatch,
305  const labelList& sourceNewNbrProc
306 )
307 {
308  Pout<< nl
309  << "Current coupling info:"
310  << endl;
311 
312  forAll(sourceFace, bFacei)
313  {
314  label meshFacei = mesh.nInternalFaces() + bFacei;
315 
316  Pout<< " meshFace:" << meshFacei
317  << " fc:" << mesh.faceCentres()[meshFacei]
318  << " connects to proc:" << sourceProc[bFacei]
319  << "/face:" << sourceFace[bFacei]
320  << " which will move to proc:" << sourceNewNbrProc[bFacei]
321  << endl;
322  }
323 }
324 
325 
326 // Finds (non-empty) patch that exposed internal and proc faces can be put into.
327 Foam::label Foam::fvMeshDistribute::findNonEmptyPatch() const
328 {
329  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
330 
331  label nonEmptyPatchi = -1;
332 
333  forAllReverse(patches, patchi)
334  {
335  const polyPatch& pp = patches[patchi];
336 
337  if (!isA<emptyPolyPatch>(pp) && !pp.coupled())
338  {
339  nonEmptyPatchi = patchi;
340  break;
341  }
342  }
343 
344  if (nonEmptyPatchi == -1)
345  {
347  << "Cannot find a patch which is neither of type empty nor"
348  << " coupled in patches " << patches.names() << endl
349  << "There has to be at least one such patch for"
350  << " distribution to work" << abort(FatalError);
351  }
352 
353  if (debug)
354  {
355  Pout<< "findNonEmptyPatch : using patch " << nonEmptyPatchi
356  << " name:" << patches[nonEmptyPatchi].name()
357  << " type:" << patches[nonEmptyPatchi].type()
358  << " to put exposed faces into." << endl;
359  }
360 
361 
362  // Do additional test for processor patches intermingled with non-proc
363  // patches.
364  label procPatchi = -1;
365 
366  forAll(patches, patchi)
367  {
368  if (isA<processorPolyPatch>(patches[patchi]))
369  {
370  procPatchi = patchi;
371  }
372  else if (procPatchi != -1)
373  {
375  << "Processor patches should be at end of patch list."
376  << endl
377  << "Have processor patch " << procPatchi
378  << " followed by non-processor patch " << patchi
379  << " in patches " << patches.names()
380  << abort(FatalError);
381  }
382  }
383 
384  return nonEmptyPatchi;
385 }
386 
387 
388 // Delete all processor patches. Move any processor faces into the last
389 // non-processor patch.
390 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::deleteProcPatches
391 (
392  const label destinationPatch
393 )
394 {
395  // New patchID per boundary faces to be repatched. Is -1 (no change)
396  // or new patchID
397  labelList newPatchID(mesh_.nFaces() - mesh_.nInternalFaces(), -1);
398 
399  label nProcPatches = 0;
400 
401  forAll(mesh_.boundaryMesh(), patchi)
402  {
403  const polyPatch& pp = mesh_.boundaryMesh()[patchi];
404 
405  if (isA<processorPolyPatch>(pp))
406  {
407  if (debug)
408  {
409  Pout<< "Moving all faces of patch " << pp.name()
410  << " into patch " << destinationPatch
411  << endl;
412  }
413 
414  label offset = pp.start() - mesh_.nInternalFaces();
415 
416  forAll(pp, i)
417  {
418  newPatchID[offset+i] = destinationPatch;
419  }
420 
421  nProcPatches++;
422  }
423  }
424 
425  // Note: order of boundary faces been kept the same since the
426  // destinationPatch is at the end and we have visited the patches in
427  // incremental order.
428  labelListList dummyFaceMaps;
429  autoPtr<mapPolyMesh> map = repatch(newPatchID, dummyFaceMaps);
430 
431 
432  // Delete (now empty) processor patches.
433  {
434  labelList oldToNew(identity(mesh_.boundaryMesh().size()));
435  label newI = 0;
436  // Non processor patches first
437  forAll(mesh_.boundaryMesh(), patchi)
438  {
439  if (!isA<processorPolyPatch>(mesh_.boundaryMesh()[patchi]))
440  {
441  oldToNew[patchi] = newI++;
442  }
443  }
444  label nNonProcPatches = newI;
445 
446  // Processor patches as last
447  forAll(mesh_.boundaryMesh(), patchi)
448  {
449  if (isA<processorPolyPatch>(mesh_.boundaryMesh()[patchi]))
450  {
451  oldToNew[patchi] = newI++;
452  }
453  }
454  fvMeshTools::reorderPatches(mesh_, oldToNew, nNonProcPatches, false);
455  }
456 
457  return map;
458 }
459 
460 
461 // Repatch the mesh.
462 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::repatch
463 (
464  const labelList& newPatchID, // per boundary face -1 or new patchID
465  labelListList& constructFaceMap
466 )
467 {
468  polyTopoChange meshMod(mesh_);
469 
470  forAll(newPatchID, bFacei)
471  {
472  if (newPatchID[bFacei] != -1)
473  {
474  label facei = mesh_.nInternalFaces() + bFacei;
475 
476  label zoneID = mesh_.faceZones().whichZone(facei);
477  bool zoneFlip = false;
478 
479  if (zoneID >= 0)
480  {
481  const faceZone& fZone = mesh_.faceZones()[zoneID];
482  zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
483  }
484 
485  meshMod.setAction
486  (
488  (
489  mesh_.faces()[facei], // modified face
490  facei, // label of face
491  mesh_.faceOwner()[facei], // owner
492  -1, // neighbour
493  false, // face flip
494  newPatchID[bFacei], // patch for face
495  false, // remove from zone
496  zoneID, // zone for face
497  zoneFlip // face flip in zone
498  )
499  );
500  }
501  }
502 
503 
504  // Do mapping of fields from one patchField to the other ourselves since
505  // is currently not supported by updateMesh.
506 
507  // Store boundary fields (we only do this for surfaceFields)
509  saveBoundaryFields<scalar, surfaceMesh>(sFlds);
511  saveBoundaryFields<vector, surfaceMesh>(vFlds);
513  saveBoundaryFields<sphericalTensor, surfaceMesh>(sptFlds);
515  saveBoundaryFields<symmTensor, surfaceMesh>(sytFlds);
517  saveBoundaryFields<tensor, surfaceMesh>(tFlds);
518 
519  // Change the mesh (no inflation). Note: parallel comms allowed.
520  //
521  // NOTE: there is one very particular problem with this ordering.
522  // We first create the processor patches and use these to merge out
523  // shared points (see mergeSharedPoints below). So temporarily points
524  // and edges do not match!
525 
526  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
527 
528  // Update fields. No inflation, parallel sync.
529  mesh_.updateMesh(map);
530 
531  // Map patch fields using stored boundary fields. Note: assumes order
532  // of fields has not changed in object registry!
533  mapBoundaryFields<scalar, surfaceMesh>(map, sFlds);
534  mapBoundaryFields<vector, surfaceMesh>(map, vFlds);
535  mapBoundaryFields<sphericalTensor, surfaceMesh>(map, sptFlds);
536  mapBoundaryFields<symmTensor, surfaceMesh>(map, sytFlds);
537  mapBoundaryFields<tensor, surfaceMesh>(map, tFlds);
538 
539 
540  // Move mesh (since morphing does not do this)
541  if (map().hasMotionPoints())
542  {
543  mesh_.movePoints(map().preMotionPoints());
544  }
545 
546  // Adapt constructMaps.
547 
548  if (debug)
549  {
550  label index = findIndex(map().reverseFaceMap(), -1);
551 
552  if (index != -1)
553  {
555  << "reverseFaceMap contains -1 at index:"
556  << index << endl
557  << "This means that the repatch operation was not just"
558  << " a shuffle?" << abort(FatalError);
559  }
560  }
561 
562  forAll(constructFaceMap, proci)
563  {
564  inplaceRenumberWithFlip
565  (
566  map().reverseFaceMap(),
567  false,
568  true,
569  constructFaceMap[proci]
570  );
571  }
572 
573 
574  return map;
575 }
576 
577 
578 // Detect shared points. Need processor patches to be present.
579 // Background: when adding bits of mesh one can get points which
580 // share the same position but are only detectable to be topologically
581 // the same point when doing parallel analysis. This routine will
582 // merge those points.
583 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
584 (
585  labelListList& constructPointMap
586 )
587 {
588  // Find out which sets of points get merged and create a map from
589  // mesh point to unique point.
590  Map<label> pointToMaster
591  (
593  (
594  mesh_,
595  mergeTol_
596  )
597  );
598 
599  if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
600  {
601  return autoPtr<mapPolyMesh>(NULL);
602  }
603 
604  polyTopoChange meshMod(mesh_);
605 
606  fvMeshAdder::mergePoints(mesh_, pointToMaster, meshMod);
607 
608  // Change the mesh (no inflation). Note: parallel comms allowed.
609  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
610 
611  // Update fields. No inflation, parallel sync.
612  mesh_.updateMesh(map);
613 
614  // Adapt constructMaps for merged points.
615  forAll(constructPointMap, proci)
616  {
617  labelList& constructMap = constructPointMap[proci];
618 
619  forAll(constructMap, i)
620  {
621  label oldPointi = constructMap[i];
622 
623  label newPointi = map().reversePointMap()[oldPointi];
624 
625  if (newPointi < -1)
626  {
627  constructMap[i] = -newPointi-2;
628  }
629  else if (newPointi >= 0)
630  {
631  constructMap[i] = newPointi;
632  }
633  else
634  {
636  << "Problem. oldPointi:" << oldPointi
637  << " newPointi:" << newPointi << abort(FatalError);
638  }
639  }
640  }
641  return map;
642 }
643 
644 
645 // Construct the local environment of all boundary faces.
646 void Foam::fvMeshDistribute::getNeighbourData
647 (
648  const labelList& distribution,
649  labelList& sourceFace,
650  labelList& sourceProc,
651  labelList& sourcePatch,
652  labelList& sourceNewNbrProc
653 ) const
654 {
655  label nBnd = mesh_.nFaces() - mesh_.nInternalFaces();
656  sourceFace.setSize(nBnd);
657  sourceProc.setSize(nBnd);
658  sourcePatch.setSize(nBnd);
659  sourceNewNbrProc.setSize(nBnd);
660 
661  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
662 
663  // Get neighbouring meshFace labels and new processor of coupled boundaries.
664  labelList nbrFaces(nBnd, -1);
665  labelList nbrNewNbrProc(nBnd, -1);
666 
667  forAll(patches, patchi)
668  {
669  const polyPatch& pp = patches[patchi];
670 
671  if (pp.coupled())
672  {
673  label offset = pp.start() - mesh_.nInternalFaces();
674 
675  // Mesh labels of faces on this side
676  forAll(pp, i)
677  {
678  label bndI = offset + i;
679  nbrFaces[bndI] = pp.start()+i;
680  }
681 
682  // Which processor they will end up on
683  SubList<label>(nbrNewNbrProc, pp.size(), offset) =
684  UIndirectList<label>(distribution, pp.faceCells())();
685  }
686  }
687 
688 
689  // Exchange the boundary data
690  syncTools::swapBoundaryFaceList(mesh_, nbrFaces);
691  syncTools::swapBoundaryFaceList(mesh_, nbrNewNbrProc);
692 
693 
694  forAll(patches, patchi)
695  {
696  const polyPatch& pp = patches[patchi];
697  label offset = pp.start() - mesh_.nInternalFaces();
698 
699  if (isA<processorPolyPatch>(pp))
700  {
701  const processorPolyPatch& procPatch =
702  refCast<const processorPolyPatch>(pp);
703 
704  // Check which of the two faces we store.
705 
706  if (procPatch.owner())
707  {
708  // Use my local face labels
709  forAll(pp, i)
710  {
711  label bndI = offset + i;
712  sourceFace[bndI] = pp.start()+i;
713  sourceProc[bndI] = Pstream::myProcNo();
714  sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
715  }
716  }
717  else
718  {
719  // Use my neighbours face labels
720  forAll(pp, i)
721  {
722  label bndI = offset + i;
723  sourceFace[bndI] = nbrFaces[bndI];
724  sourceProc[bndI] = procPatch.neighbProcNo();
725  sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
726  }
727  }
728 
729 
730  label patchi = -1;
731  if (isA<processorCyclicPolyPatch>(pp))
732  {
733  patchi = refCast<const processorCyclicPolyPatch>
734  (
735  pp
736  ).referPatchID();
737  }
738 
739  forAll(pp, i)
740  {
741  label bndI = offset + i;
742  sourcePatch[bndI] = patchi;
743  }
744  }
745  else if (isA<cyclicPolyPatch>(pp))
746  {
747  const cyclicPolyPatch& cpp = refCast<const cyclicPolyPatch>(pp);
748 
749  if (cpp.owner())
750  {
751  forAll(pp, i)
752  {
753  label bndI = offset + i;
754  sourceFace[bndI] = pp.start()+i;
755  sourceProc[bndI] = Pstream::myProcNo();
756  sourcePatch[bndI] = patchi;
757  sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
758  }
759  }
760  else
761  {
762  forAll(pp, i)
763  {
764  label bndI = offset + i;
765  sourceFace[bndI] = nbrFaces[bndI];
766  sourceProc[bndI] = Pstream::myProcNo();
767  sourcePatch[bndI] = patchi;
768  sourceNewNbrProc[bndI] = nbrNewNbrProc[bndI];
769  }
770  }
771  }
772  else
773  {
774  // Normal (physical) boundary
775  forAll(pp, i)
776  {
777  label bndI = offset + i;
778  sourceFace[bndI] = -1;
779  sourceProc[bndI] = -1;
780  sourcePatch[bndI] = patchi;
781  sourceNewNbrProc[bndI] = -1;
782  }
783  }
784  }
785 }
786 
787 
788 // Subset the neighbourCell/neighbourProc fields
789 void Foam::fvMeshDistribute::subsetBoundaryData
790 (
791  const fvMesh& mesh,
792  const labelList& faceMap,
793  const labelList& cellMap,
794 
795  const labelList& oldDistribution,
796  const labelList& oldFaceOwner,
797  const labelList& oldFaceNeighbour,
798  const label oldInternalFaces,
799 
800  const labelList& sourceFace,
801  const labelList& sourceProc,
802  const labelList& sourcePatch,
803  const labelList& sourceNewNbrProc,
804 
805  labelList& subFace,
806  labelList& subProc,
807  labelList& subPatch,
808  labelList& subNewNbrProc
809 )
810 {
811  subFace.setSize(mesh.nFaces() - mesh.nInternalFaces());
812  subProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
813  subPatch.setSize(mesh.nFaces() - mesh.nInternalFaces());
814  subNewNbrProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
815 
816  forAll(subFace, newBFacei)
817  {
818  label newFacei = newBFacei + mesh.nInternalFaces();
819 
820  label oldFacei = faceMap[newFacei];
821 
822  // Was oldFacei internal face? If so which side did we get.
823  if (oldFacei < oldInternalFaces)
824  {
825  subFace[newBFacei] = oldFacei;
826  subProc[newBFacei] = Pstream::myProcNo();
827  subPatch[newBFacei] = -1;
828 
829  label oldOwn = oldFaceOwner[oldFacei];
830  label oldNei = oldFaceNeighbour[oldFacei];
831 
832  if (oldOwn == cellMap[mesh.faceOwner()[newFacei]])
833  {
834  // We kept the owner side. Where does the neighbour move to?
835  subNewNbrProc[newBFacei] = oldDistribution[oldNei];
836  }
837  else
838  {
839  // We kept the neighbour side.
840  subNewNbrProc[newBFacei] = oldDistribution[oldOwn];
841  }
842  }
843  else
844  {
845  // Was boundary face. Take over boundary information
846  label oldBFacei = oldFacei - oldInternalFaces;
847 
848  subFace[newBFacei] = sourceFace[oldBFacei];
849  subProc[newBFacei] = sourceProc[oldBFacei];
850  subPatch[newBFacei] = sourcePatch[oldBFacei];
851  subNewNbrProc[newBFacei] = sourceNewNbrProc[oldBFacei];
852  }
853  }
854 }
855 
856 
857 // Find cells on mesh whose faceID/procID match the neighbour cell/proc of
858 // domainMesh. Store the matching face.
859 void Foam::fvMeshDistribute::findCouples
860 (
861  const primitiveMesh& mesh,
862  const labelList& sourceFace,
863  const labelList& sourceProc,
864  const labelList& sourcePatch,
865 
866  const label domain,
867  const primitiveMesh& domainMesh,
868  const labelList& domainFace,
869  const labelList& domainProc,
870  const labelList& domainPatch,
871 
872  labelList& masterCoupledFaces,
873  labelList& slaveCoupledFaces
874 )
875 {
876  // Store domain neighbour as map so we can easily look for pair
877  // with same face+proc.
879 
880  forAll(domainProc, bFacei)
881  {
882  if (domainProc[bFacei] != -1 && domainPatch[bFacei] == -1)
883  {
884  map.insert
885  (
886  labelPair(domainFace[bFacei], domainProc[bFacei]),
887  bFacei
888  );
889  }
890  }
891 
892 
893  // Try to match mesh data.
894 
895  masterCoupledFaces.setSize(domainFace.size());
896  slaveCoupledFaces.setSize(domainFace.size());
897  label coupledI = 0;
898 
899  forAll(sourceFace, bFacei)
900  {
901  if (sourceProc[bFacei] != -1 && sourcePatch[bFacei] == -1)
902  {
903  labelPair myData(sourceFace[bFacei], sourceProc[bFacei]);
904 
906  iter = map.find(myData);
907 
908  if (iter != map.end())
909  {
910  label nbrBFacei = iter();
911 
912  masterCoupledFaces[coupledI] = mesh.nInternalFaces() + bFacei;
913  slaveCoupledFaces[coupledI] =
914  domainMesh.nInternalFaces()
915  + nbrBFacei;
916 
917  coupledI++;
918  }
919  }
920  }
921 
922  masterCoupledFaces.setSize(coupledI);
923  slaveCoupledFaces.setSize(coupledI);
924 
925  if (debug)
926  {
927  Pout<< "findCouples : found " << coupledI
928  << " faces that will be stitched" << nl << endl;
929  }
930 }
931 
932 
933 // Map data on boundary faces to new mesh (resulting from adding two meshes)
934 Foam::labelList Foam::fvMeshDistribute::mapBoundaryData
935 (
936  const primitiveMesh& mesh, // mesh after adding
937  const mapAddedPolyMesh& map,
938  const labelList& boundaryData0, // mesh before adding
939  const label nInternalFaces1,
940  const labelList& boundaryData1 // added mesh
941 )
942 {
943  labelList newBoundaryData(mesh.nFaces() - mesh.nInternalFaces());
944 
945  forAll(boundaryData0, oldBFacei)
946  {
947  label newFacei = map.oldFaceMap()[oldBFacei + map.nOldInternalFaces()];
948 
949  // Face still exists (is necessary?) and still boundary face
950  if (newFacei >= 0 && newFacei >= mesh.nInternalFaces())
951  {
952  newBoundaryData[newFacei - mesh.nInternalFaces()] =
953  boundaryData0[oldBFacei];
954  }
955  }
956 
957  forAll(boundaryData1, addedBFacei)
958  {
959  label newFacei = map.addedFaceMap()[addedBFacei + nInternalFaces1];
960 
961  if (newFacei >= 0 && newFacei >= mesh.nInternalFaces())
962  {
963  newBoundaryData[newFacei - mesh.nInternalFaces()] =
964  boundaryData1[addedBFacei];
965  }
966  }
967 
968  return newBoundaryData;
969 }
970 
971 
972 // Remove cells. Add all exposed faces to patch oldInternalPatchi
973 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::doRemoveCells
974 (
975  const labelList& cellsToRemove,
976  const label oldInternalPatchi
977 )
978 {
979  // Mesh change engine
980  polyTopoChange meshMod(mesh_);
981 
982  // Cell removal topo engine. Do NOT synchronize parallel since
983  // we are doing a local cell removal.
984  removeCells cellRemover(mesh_, false);
985 
986  // Get all exposed faces
987  labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
988 
989  // Insert the topo changes
990  cellRemover.setRefinement
991  (
992  cellsToRemove,
993  exposedFaces,
994  labelList(exposedFaces.size(), oldInternalPatchi), // patch for exposed
995  // faces.
996  meshMod
997  );
998 
999 
1001  //tmp<surfaceScalarField> sfld(generateTestField(mesh_));
1002 
1003  // Save internal fields (note: not as DimensionedFields since would
1004  // get mapped)
1005  PtrList<Field<scalar>> sFlds;
1006  saveInternalFields(sFlds);
1007  PtrList<Field<vector>> vFlds;
1008  saveInternalFields(vFlds);
1010  saveInternalFields(sptFlds);
1011  PtrList<Field<symmTensor>> sytFlds;
1012  saveInternalFields(sytFlds);
1013  PtrList<Field<tensor>> tFlds;
1014  saveInternalFields(tFlds);
1015 
1016  // Change the mesh. No inflation. Note: no parallel comms allowed.
1017  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, false);
1018 
1019  // Update fields
1020  mesh_.updateMesh(map);
1021 
1022 
1023  // Any exposed faces in a surfaceField will not be mapped. Map the value
1024  // of these separately (until there is support in all PatchFields for
1025  // mapping from internal faces ...)
1026 
1027  mapExposedFaces(map(), sFlds);
1028  mapExposedFaces(map(), vFlds);
1029  mapExposedFaces(map(), sptFlds);
1030  mapExposedFaces(map(), sytFlds);
1031  mapExposedFaces(map(), tFlds);
1032 
1033 
1035  //testField(sfld);
1036 
1037 
1038  // Move mesh (since morphing does not do this)
1039  if (map().hasMotionPoints())
1040  {
1041  mesh_.movePoints(map().preMotionPoints());
1042  }
1043 
1044  return map;
1045 }
1046 
1047 
1048 // Delete and add processor patches. Changes mesh and returns per neighbour proc
1049 // the processor patchID.
1050 void Foam::fvMeshDistribute::addProcPatches
1051 (
1052  const labelList& nbrProc, // processor that neighbour is now on
1053  const labelList& referPatchID, // patchID (or -1) I originated from
1054  List<Map<label>>& procPatchID
1055 )
1056 {
1057  // Now use the neighbourFace/Proc to repatch the mesh. These lists
1058  // contain for all current boundary faces the global patchID (for non-proc
1059  // patch) or the processor.
1060 
1061  // Determine a visit order such that the processor patches get added
1062  // in order of increasing neighbour processor (and for same neighbour
1063  // processor (in case of processor cyclics) in order of increasing
1064  // 'refer' patch)
1065  labelList indices;
1066  sortedOrder(nbrProc, indices, lessProcPatches(nbrProc, referPatchID));
1067 
1068  procPatchID.setSize(Pstream::nProcs());
1069 
1070  forAll(indices, i)
1071  {
1072  label bFacei = indices[i];
1073  label proci = nbrProc[bFacei];
1074 
1075  if (proci != -1 && proci != Pstream::myProcNo())
1076  {
1077  if (!procPatchID[proci].found(referPatchID[bFacei]))
1078  {
1079  // No patch for neighbour yet. Is either a normal processor
1080  // patch or a processorCyclic patch.
1081 
1082  if (referPatchID[bFacei] == -1)
1083  {
1084  // Ordinary processor boundary
1085 
1087  (
1088  0, // size
1089  mesh_.nFaces(),
1090  mesh_.boundaryMesh().size(),
1091  mesh_.boundaryMesh(),
1093  proci
1094  );
1095 
1096  procPatchID[proci].insert
1097  (
1098  referPatchID[bFacei],
1100  (
1101  mesh_,
1102  pp,
1103  dictionary(), // optional per field patchField
1105  false // not parallel sync
1106  )
1107  );
1108  }
1109  else
1110  {
1111  const coupledPolyPatch& pcPatch
1112  = refCast<const coupledPolyPatch>
1113  (
1114  mesh_.boundaryMesh()[referPatchID[bFacei]]
1115  );
1117  (
1118  0, // size
1119  mesh_.nFaces(),
1120  mesh_.boundaryMesh().size(),
1121  mesh_.boundaryMesh(),
1123  proci,
1124  pcPatch.name(),
1125  pcPatch.transform()
1126  );
1127 
1128  procPatchID[proci].insert
1129  (
1130  referPatchID[bFacei],
1132  (
1133  mesh_,
1134  pp,
1135  dictionary(), // optional per field patchField
1137  false // not parallel sync
1138  )
1139  );
1140  }
1141  }
1142  }
1143  }
1144 }
1145 
1146 
1147 // Get boundary faces to be repatched. Is -1 or new patchID
1148 Foam::labelList Foam::fvMeshDistribute::getBoundaryPatch
1149 (
1150  const labelList& nbrProc, // new processor per boundary face
1151  const labelList& referPatchID, // patchID (or -1) I originated from
1152  const List<Map<label>>& procPatchID // per proc the new procPatches
1153 )
1154 {
1155  labelList patchIDs(nbrProc);
1156 
1157  forAll(nbrProc, bFacei)
1158  {
1159  if (nbrProc[bFacei] == Pstream::myProcNo())
1160  {
1161  label origPatchi = referPatchID[bFacei];
1162  patchIDs[bFacei] = origPatchi;
1163  }
1164  else if (nbrProc[bFacei] != -1)
1165  {
1166  label origPatchi = referPatchID[bFacei];
1167  patchIDs[bFacei] = procPatchID[nbrProc[bFacei]][origPatchi];
1168  }
1169  else
1170  {
1171  patchIDs[bFacei] = -1;
1172  }
1173  }
1174  return patchIDs;
1175 }
1176 
1177 
1178 // Send mesh and coupling data.
1179 void Foam::fvMeshDistribute::sendMesh
1180 (
1181  const label domain,
1182  const fvMesh& mesh,
1183 
1184  const wordList& pointZoneNames,
1185  const wordList& faceZoneNames,
1186  const wordList& cellZoneNames,
1187 
1188  const labelList& sourceFace,
1189  const labelList& sourceProc,
1190  const labelList& sourcePatch,
1191  const labelList& sourceNewNbrProc,
1192  Ostream& toDomain
1193 )
1194 {
1195  if (debug)
1196  {
1197  Pout<< "Sending to domain " << domain << nl
1198  << " nPoints:" << mesh.nPoints() << nl
1199  << " nFaces:" << mesh.nFaces() << nl
1200  << " nCells:" << mesh.nCells() << nl
1201  << " nPatches:" << mesh.boundaryMesh().size() << nl
1202  << endl;
1203  }
1204 
1205  // Assume sparse, possibly overlapping point zones. Get contents
1206  // in merged-zone indices.
1207  CompactListList<label> zonePoints;
1208  {
1209  const pointZoneMesh& pointZones = mesh.pointZones();
1210 
1211  labelList rowSizes(pointZoneNames.size(), 0);
1212 
1213  forAll(pointZoneNames, nameI)
1214  {
1215  label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
1216 
1217  if (myZoneID != -1)
1218  {
1219  rowSizes[nameI] = pointZones[myZoneID].size();
1220  }
1221  }
1222  zonePoints.setSize(rowSizes);
1223 
1224  forAll(pointZoneNames, nameI)
1225  {
1226  label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
1227 
1228  if (myZoneID != -1)
1229  {
1230  zonePoints[nameI].deepCopy(pointZones[myZoneID]);
1231  }
1232  }
1233  }
1234 
1235  // Assume sparse, possibly overlapping face zones
1236  CompactListList<label> zoneFaces;
1237  CompactListList<bool> zoneFaceFlip;
1238  {
1239  const faceZoneMesh& faceZones = mesh.faceZones();
1240 
1241  labelList rowSizes(faceZoneNames.size(), 0);
1242 
1243  forAll(faceZoneNames, nameI)
1244  {
1245  label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
1246 
1247  if (myZoneID != -1)
1248  {
1249  rowSizes[nameI] = faceZones[myZoneID].size();
1250  }
1251  }
1252 
1253  zoneFaces.setSize(rowSizes);
1254  zoneFaceFlip.setSize(rowSizes);
1255 
1256  forAll(faceZoneNames, nameI)
1257  {
1258  label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
1259 
1260  if (myZoneID != -1)
1261  {
1262  zoneFaces[nameI].deepCopy(faceZones[myZoneID]);
1263  zoneFaceFlip[nameI].deepCopy(faceZones[myZoneID].flipMap());
1264  }
1265  }
1266  }
1267 
1268  // Assume sparse, possibly overlapping cell zones
1269  CompactListList<label> zoneCells;
1270  {
1271  const cellZoneMesh& cellZones = mesh.cellZones();
1272 
1273  labelList rowSizes(cellZoneNames.size(), 0);
1274 
1275  forAll(cellZoneNames, nameI)
1276  {
1277  label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
1278 
1279  if (myZoneID != -1)
1280  {
1281  rowSizes[nameI] = cellZones[myZoneID].size();
1282  }
1283  }
1284 
1285  zoneCells.setSize(rowSizes);
1286 
1287  forAll(cellZoneNames, nameI)
1288  {
1289  label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
1290 
1291  if (myZoneID != -1)
1292  {
1293  zoneCells[nameI].deepCopy(cellZones[myZoneID]);
1294  }
1295  }
1296  }
1298  //labelList cellZoneID;
1299  //if (hasCellZones)
1300  //{
1301  // cellZoneID.setSize(mesh.nCells());
1302  // cellZoneID = -1;
1303  //
1304  // const cellZoneMesh& cellZones = mesh.cellZones();
1305  //
1306  // forAll(cellZones, zoneI)
1307  // {
1308  // UIndirectList<label>(cellZoneID, cellZones[zoneI]) = zoneI;
1309  // }
1310  //}
1311 
1312  // Send
1313  toDomain
1314  << mesh.points()
1316  << mesh.faceOwner()
1317  << mesh.faceNeighbour()
1318  << mesh.boundaryMesh()
1319 
1320  << zonePoints
1321  << zoneFaces
1322  << zoneFaceFlip
1323  << zoneCells
1324 
1325  << sourceFace
1326  << sourceProc
1327  << sourcePatch
1328  << sourceNewNbrProc;
1329 
1330 
1331  if (debug)
1332  {
1333  Pout<< "Started sending mesh to domain " << domain
1334  << endl;
1335  }
1336 }
1337 
1338 
1339 // Receive mesh. Opposite of sendMesh
1340 Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
1341 (
1342  const label domain,
1343  const wordList& pointZoneNames,
1344  const wordList& faceZoneNames,
1345  const wordList& cellZoneNames,
1346  const Time& runTime,
1347  labelList& domainSourceFace,
1348  labelList& domainSourceProc,
1349  labelList& domainSourcePatch,
1350  labelList& domainSourceNewNbrProc,
1351  Istream& fromNbr
1352 )
1353 {
1354  pointField domainPoints(fromNbr);
1355  faceList domainFaces = CompactListList<label, face>(fromNbr)();
1356  labelList domainAllOwner(fromNbr);
1357  labelList domainAllNeighbour(fromNbr);
1358  PtrList<entry> patchEntries(fromNbr);
1359 
1360  CompactListList<label> zonePoints(fromNbr);
1361  CompactListList<label> zoneFaces(fromNbr);
1362  CompactListList<bool> zoneFaceFlip(fromNbr);
1363  CompactListList<label> zoneCells(fromNbr);
1364 
1365  fromNbr
1366  >> domainSourceFace
1367  >> domainSourceProc
1368  >> domainSourcePatch
1369  >> domainSourceNewNbrProc;
1370 
1371  // Construct fvMesh
1372  autoPtr<fvMesh> domainMeshPtr
1373  (
1374  new fvMesh
1375  (
1376  IOobject
1377  (
1379  runTime.timeName(),
1380  runTime,
1382  ),
1383  xferMove(domainPoints),
1384  xferMove(domainFaces),
1385  xferMove(domainAllOwner),
1386  xferMove(domainAllNeighbour),
1387  false // no parallel comms
1388  )
1389  );
1390  fvMesh& domainMesh = domainMeshPtr();
1391 
1392  List<polyPatch*> patches(patchEntries.size());
1393 
1394  forAll(patchEntries, patchi)
1395  {
1397  (
1398  patchEntries[patchi].keyword(),
1399  patchEntries[patchi].dict(),
1400  patchi,
1401  domainMesh.boundaryMesh()
1402  ).ptr();
1403  }
1404  // Add patches; no parallel comms
1405  domainMesh.addFvPatches(patches, false);
1406 
1407  // Construct zones
1408  List<pointZone*> pZonePtrs(pointZoneNames.size());
1409  forAll(pZonePtrs, i)
1410  {
1411  pZonePtrs[i] = new pointZone
1412  (
1413  pointZoneNames[i],
1414  zonePoints[i],
1415  i,
1416  domainMesh.pointZones()
1417  );
1418  }
1419 
1420  List<faceZone*> fZonePtrs(faceZoneNames.size());
1421  forAll(fZonePtrs, i)
1422  {
1423  fZonePtrs[i] = new faceZone
1424  (
1425  faceZoneNames[i],
1426  zoneFaces[i],
1427  zoneFaceFlip[i],
1428  i,
1429  domainMesh.faceZones()
1430  );
1431  }
1432 
1433  List<cellZone*> cZonePtrs(cellZoneNames.size());
1434  forAll(cZonePtrs, i)
1435  {
1436  cZonePtrs[i] = new cellZone
1437  (
1438  cellZoneNames[i],
1439  zoneCells[i],
1440  i,
1441  domainMesh.cellZones()
1442  );
1443  }
1444  domainMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
1445 
1446  return domainMeshPtr;
1447 }
1448 
1449 
1450 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1451 
1452 // Construct from components
1453 Foam::fvMeshDistribute::fvMeshDistribute(fvMesh& mesh, const scalar mergeTol)
1454 :
1455  mesh_(mesh),
1456  mergeTol_(mergeTol)
1457 {}
1458 
1459 
1460 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1461 
1464  const labelList& distribution
1465 )
1466 {
1467  labelList nCells(Pstream::nProcs(), 0);
1468  forAll(distribution, celli)
1469  {
1470  label newProc = distribution[celli];
1471 
1472  if (newProc < 0 || newProc >= Pstream::nProcs())
1473  {
1475  << "Distribution should be in range 0.." << Pstream::nProcs()-1
1476  << endl
1477  << "At index " << celli << " distribution:" << newProc
1478  << abort(FatalError);
1479  }
1480  nCells[newProc]++;
1481  }
1482  return nCells;
1483 }
1484 
1485 
1488  const labelList& distribution
1489 )
1490 {
1491  // Some checks on distribution
1492  if (distribution.size() != mesh_.nCells())
1493  {
1495  << "Size of distribution:"
1496  << distribution.size() << " mesh nCells:" << mesh_.nCells()
1497  << abort(FatalError);
1498  }
1499 
1500 
1501  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1502 
1503  // Check all processors have same non-proc patches in same order.
1504  if (patches.checkParallelSync(true))
1505  {
1507  << "This application requires all non-processor patches"
1508  << " to be present in the same order on all patches" << nl
1509  << "followed by the processor patches (which of course are unique)."
1510  << nl
1511  << "Local patches:" << mesh_.boundaryMesh().names()
1512  << abort(FatalError);
1513  }
1514 
1515  // Save some data for mapping later on
1516  const label nOldPoints(mesh_.nPoints());
1517  const label nOldFaces(mesh_.nFaces());
1518  const label nOldCells(mesh_.nCells());
1519  labelList oldPatchStarts(patches.size());
1520  labelList oldPatchNMeshPoints(patches.size());
1521  forAll(patches, patchi)
1522  {
1523  oldPatchStarts[patchi] = patches[patchi].start();
1524  oldPatchNMeshPoints[patchi] = patches[patchi].nPoints();
1525  }
1526 
1527 
1528 
1529  // Short circuit trivial case.
1530  if (!Pstream::parRun())
1531  {
1532  // Collect all maps and return
1534  (
1536  (
1537  mesh_,
1538 
1539  nOldPoints,
1540  nOldFaces,
1541  nOldCells,
1542  oldPatchStarts.xfer(),
1543  oldPatchNMeshPoints.xfer(),
1544 
1545  labelListList(1, identity(mesh_.nPoints())).xfer(),//subPointMap
1546  labelListList(1, identity(mesh_.nFaces())).xfer(), //subFaceMap
1547  labelListList(1, identity(mesh_.nCells())).xfer(), //subCellMap
1548  labelListList(1, identity(patches.size())).xfer(), //subPatchMap
1549 
1550  labelListList(1, identity(mesh_.nPoints())).xfer(),//pointMap
1551  labelListList(1, identity(mesh_.nFaces())).xfer(), //faceMap
1552  labelListList(1, identity(mesh_.nCells())).xfer(), //cellMap
1553  labelListList(1, identity(patches.size())).xfer() //patchMap
1554  )
1555  );
1556  }
1557 
1558 
1559  // Collect any zone names
1560  const wordList pointZoneNames(mergeWordList(mesh_.pointZones().names()));
1561  const wordList faceZoneNames(mergeWordList(mesh_.faceZones().names()));
1562  const wordList cellZoneNames(mergeWordList(mesh_.cellZones().names()));
1563 
1564 
1565 
1566  // Local environment of all boundary faces
1567  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1568 
1569  // A face is uniquely defined by
1570  // - proc
1571  // - local face no
1572  //
1573  // To glue the parts of meshes which can get sent from anywhere we
1574  // need to know on boundary faces what the above tuple on both sides is.
1575  // So we need to maintain:
1576  // - original face
1577  // - original processor id (= trivial)
1578  // For coupled boundaries (where the faces are 'duplicate') we take the
1579  // lowest numbered processor as the data to store.
1580  //
1581  // Additionally to create the procboundaries we need to know where the owner
1582  // cell on the other side now is: newNeighbourProc.
1583  //
1584 
1585  // physical boundary:
1586  // sourceProc = -1
1587  // sourceNewNbrProc = -1
1588  // sourceFace = -1
1589  // sourcePatch = patchID
1590  // processor boundary:
1591  // sourceProc = proc (on owner side)
1592  // sourceNewNbrProc = distribution of coupled cell
1593  // sourceFace = face (on owner side)
1594  // sourcePatch = -1
1595  // ?cyclic:
1596  // ? sourceProc = proc
1597  // ? sourceNewNbrProc = distribution of coupled cell
1598  // ? sourceFace = face (on owner side)
1599  // ? sourcePatch = patchID
1600  // processor-cyclic boundary:
1601  // sourceProc = proc (on owner side)
1602  // sourceNewNbrProc = distribution of coupled cell
1603  // sourceFace = face (on owner side)
1604  // sourcePatch = patchID
1605 
1606  labelList sourcePatch;
1607  labelList sourceFace;
1608  labelList sourceProc;
1609  labelList sourceNewNbrProc;
1610  getNeighbourData
1611  (
1612  distribution,
1613  sourceFace,
1614  sourceProc,
1615  sourcePatch,
1616  sourceNewNbrProc
1617  );
1618 
1619 
1620  // Remove meshPhi. Since this would otherwise disappear anyway
1621  // during topo changes and we have to guarantee that all the fields
1622  // can be sent.
1623  mesh_.clearOut();
1624  mesh_.resetMotion();
1625 
1626  // Get data to send. Make sure is synchronised
1627  const wordList volScalars(mesh_.names(volScalarField::typeName));
1628  checkEqualWordList("volScalarFields", volScalars);
1629  const wordList volVectors(mesh_.names(volVectorField::typeName));
1630  checkEqualWordList("volVectorFields", volVectors);
1631  const wordList volSphereTensors
1632  (
1634  );
1635  checkEqualWordList("volSphericalTensorFields", volSphereTensors);
1636  const wordList volSymmTensors(mesh_.names(volSymmTensorField::typeName));
1637  checkEqualWordList("volSymmTensorFields", volSymmTensors);
1638  const wordList volTensors(mesh_.names(volTensorField::typeName));
1639  checkEqualWordList("volTensorField", volTensors);
1640 
1641  const wordList surfScalars(mesh_.names(surfaceScalarField::typeName));
1642  checkEqualWordList("surfaceScalarFields", surfScalars);
1643  const wordList surfVectors(mesh_.names(surfaceVectorField::typeName));
1644  checkEqualWordList("surfaceVectorFields", surfVectors);
1645  const wordList surfSphereTensors
1646  (
1648  );
1649  checkEqualWordList("surfaceSphericalTensorFields", surfSphereTensors);
1650  const wordList surfSymmTensors
1651  (
1653  );
1654  checkEqualWordList("surfaceSymmTensorFields", surfSymmTensors);
1655  const wordList surfTensors(mesh_.names(surfaceTensorField::typeName));
1656  checkEqualWordList("surfaceTensorFields", surfTensors);
1657 
1658  typedef volScalarField::Internal dimScalType;
1659  const wordList dimScalars(mesh_.names(dimScalType::typeName));
1660  checkEqualWordList("volScalarField::Internal", dimScalars);
1661 
1662  typedef volVectorField::Internal dimVecType;
1663  const wordList dimVectors(mesh_.names(dimVecType::typeName));
1664  checkEqualWordList("volVectorField::Internal", dimVectors);
1665 
1666  typedef volSphericalTensorField::Internal dimSphereType;
1667  const wordList dimSphereTensors(mesh_.names(dimSphereType::typeName));
1668  checkEqualWordList
1669  (
1670  "volSphericalTensorField::Internal",
1671  dimSphereTensors
1672  );
1673 
1674  typedef volSymmTensorField::Internal dimSymmTensorType;
1675  const wordList dimSymmTensors(mesh_.names(dimSymmTensorType::typeName));
1676  checkEqualWordList
1677  (
1678  "volSymmTensorField::Internal",
1679  dimSymmTensors
1680  );
1681 
1682  typedef volTensorField::Internal dimTensorType;
1683  const wordList dimTensors(mesh_.names(dimTensorType::typeName));
1684  checkEqualWordList("volTensorField::Internal", dimTensors);
1685 
1686 
1687 
1688  // Find patch to temporarily put exposed and processor faces into.
1689  label oldInternalPatchi = findNonEmptyPatch();
1690 
1691 
1692 
1693  // Delete processor patches, starting from the back. Move all faces into
1694  // oldInternalPatchi.
1695  labelList repatchFaceMap;
1696  {
1697  autoPtr<mapPolyMesh> repatchMap = deleteProcPatches(oldInternalPatchi);
1698 
1699  // Store face map (only face ordering that changed)
1700  repatchFaceMap = repatchMap().faceMap();
1701 
1702  // Reorder all boundary face data (sourceProc, sourceFace etc.)
1703  labelList bFaceMap
1704  (
1706  (
1707  repatchMap().reverseFaceMap(),
1708  mesh_.nFaces() - mesh_.nInternalFaces(),
1709  mesh_.nInternalFaces()
1710  )
1711  - mesh_.nInternalFaces()
1712  );
1713 
1714  inplaceReorder(bFaceMap, sourceFace);
1715  inplaceReorder(bFaceMap, sourceProc);
1716  inplaceReorder(bFaceMap, sourcePatch);
1717  inplaceReorder(bFaceMap, sourceNewNbrProc);
1718  }
1719 
1720 
1721 
1722  // Print a bit.
1723  if (debug)
1724  {
1725  Pout<< nl << "MESH WITH PROC PATCHES DELETED:" << endl;
1726  printMeshInfo(mesh_);
1727  printFieldInfo<volScalarField>(mesh_);
1728  printFieldInfo<volVectorField>(mesh_);
1729  printFieldInfo<volSphericalTensorField>(mesh_);
1730  printFieldInfo<volSymmTensorField>(mesh_);
1731  printFieldInfo<volTensorField>(mesh_);
1732  printFieldInfo<surfaceScalarField>(mesh_);
1733  printFieldInfo<surfaceVectorField>(mesh_);
1734  printFieldInfo<surfaceSphericalTensorField>(mesh_);
1735  printFieldInfo<surfaceSymmTensorField>(mesh_);
1736  printFieldInfo<surfaceTensorField>(mesh_);
1737  Pout<< nl << endl;
1738  }
1739 
1740 
1741 
1742  // Maps from subsetted mesh (that is sent) back to original maps
1743  labelListList subCellMap(Pstream::nProcs());
1744  labelListList subFaceMap(Pstream::nProcs());
1745  labelListList subPointMap(Pstream::nProcs());
1746  labelListList subPatchMap(Pstream::nProcs());
1747  // Maps from subsetted mesh to reconstructed mesh
1748  labelListList constructCellMap(Pstream::nProcs());
1749  labelListList constructFaceMap(Pstream::nProcs());
1750  labelListList constructPointMap(Pstream::nProcs());
1751  labelListList constructPatchMap(Pstream::nProcs());
1752 
1753 
1754 
1755 
1756  // Find out schedule
1757  // ~~~~~~~~~~~~~~~~~
1758 
1759  labelListList nSendCells(Pstream::nProcs());
1760  nSendCells[Pstream::myProcNo()] = countCells(distribution);
1761  Pstream::gatherList(nSendCells);
1762  Pstream::scatterList(nSendCells);
1763 
1764 
1765  // Allocate buffers
1767 
1768 
1769  // What to send to neighbouring domains
1770  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1771 
1772  bool oldParRun = UPstream::parRun();
1773  UPstream::parRun() = false;
1774 
1775  forAll(nSendCells[Pstream::myProcNo()], recvProc)
1776  {
1777  if
1778  (
1779  recvProc != Pstream::myProcNo()
1780  && nSendCells[Pstream::myProcNo()][recvProc] > 0
1781  )
1782  {
1783  // Send to recvProc
1784 
1785  if (debug)
1786  {
1787  Pout<< nl
1788  << "SUBSETTING FOR DOMAIN " << recvProc
1789  << " cells to send:"
1790  << nSendCells[Pstream::myProcNo()][recvProc]
1791  << nl << endl;
1792  }
1793 
1794  // Pstream for sending mesh and fields
1795  //OPstream str(Pstream::blocking, recvProc);
1796  UOPstream str(recvProc, pBufs);
1797 
1798  // Mesh subsetting engine
1799  fvMeshSubset subsetter(mesh_);
1800 
1801  // Subset the cells of the current domain.
1802  subsetter.setLargeCellSubset
1803  (
1804  distribution,
1805  recvProc,
1806  oldInternalPatchi, // oldInternalFaces patch
1807  false // no parallel sync
1808  );
1809 
1810  subCellMap[recvProc] = subsetter.cellMap();
1811  subFaceMap[recvProc] = subsetter.faceFlipMap();
1812  inplaceRenumberWithFlip
1813  (
1814  repatchFaceMap,
1815  false, // oldToNew has flip
1816  true, // subFaceMap has flip
1817  subFaceMap[recvProc]
1818  );
1819  subPointMap[recvProc] = subsetter.pointMap();
1820  subPatchMap[recvProc] = subsetter.patchMap();
1821 
1822 
1823  // Subset the boundary fields (owner/neighbour/processor)
1824  labelList procSourceFace;
1825  labelList procSourceProc;
1826  labelList procSourcePatch;
1827  labelList procSourceNewNbrProc;
1828 
1829  subsetBoundaryData
1830  (
1831  subsetter.subMesh(),
1832  subsetter.faceMap(), // from subMesh to mesh
1833  subsetter.cellMap(), // ,, ,,
1834 
1835  distribution, // old mesh distribution
1836  mesh_.faceOwner(), // old owner
1837  mesh_.faceNeighbour(),
1838  mesh_.nInternalFaces(),
1839 
1840  sourceFace,
1841  sourceProc,
1842  sourcePatch,
1843  sourceNewNbrProc,
1844 
1845  procSourceFace,
1846  procSourceProc,
1847  procSourcePatch,
1848  procSourceNewNbrProc
1849  );
1850 
1851 
1852 
1853  // Send to neighbour
1854  sendMesh
1855  (
1856  recvProc,
1857  subsetter.subMesh(),
1858 
1859  pointZoneNames,
1860  faceZoneNames,
1861  cellZoneNames,
1862 
1863  procSourceFace,
1864  procSourceProc,
1865  procSourcePatch,
1866  procSourceNewNbrProc,
1867  str
1868  );
1869 
1870  // volFields
1871  sendFields<volScalarField>(recvProc, volScalars, subsetter, str);
1872  sendFields<volVectorField>(recvProc, volVectors, subsetter, str);
1873  sendFields<volSphericalTensorField>
1874  (
1875  recvProc,
1876  volSphereTensors,
1877  subsetter,
1878  str
1879  );
1880  sendFields<volSymmTensorField>
1881  (
1882  recvProc,
1883  volSymmTensors,
1884  subsetter,
1885  str
1886  );
1887  sendFields<volTensorField>(recvProc, volTensors, subsetter, str);
1888 
1889  // surfaceFields
1890  sendFields<surfaceScalarField>
1891  (
1892  recvProc,
1893  surfScalars,
1894  subsetter,
1895  str
1896  );
1897  sendFields<surfaceVectorField>
1898  (
1899  recvProc,
1900  surfVectors,
1901  subsetter,
1902  str
1903  );
1904  sendFields<surfaceSphericalTensorField>
1905  (
1906  recvProc,
1907  surfSphereTensors,
1908  subsetter,
1909  str
1910  );
1911  sendFields<surfaceSymmTensorField>
1912  (
1913  recvProc,
1914  surfSymmTensors,
1915  subsetter,
1916  str
1917  );
1918  sendFields<surfaceTensorField>
1919  (
1920  recvProc,
1921  surfTensors,
1922  subsetter,
1923  str
1924  );
1925 
1926  // dimensionedFields
1927  sendFields<volScalarField::Internal>
1928  (
1929  recvProc,
1930  dimScalars,
1931  subsetter,
1932  str
1933  );
1934  sendFields<volVectorField::Internal>
1935  (
1936  recvProc,
1937  dimVectors,
1938  subsetter,
1939  str
1940  );
1941  sendFields<volSphericalTensorField::Internal>
1942  (
1943  recvProc,
1944  dimSphereTensors,
1945  subsetter,
1946  str
1947  );
1948  sendFields<volSymmTensorField::Internal>
1949  (
1950  recvProc,
1951  dimSymmTensors,
1952  subsetter,
1953  str
1954  );
1955  sendFields<volTensorField::Internal>
1956  (
1957  recvProc,
1958  dimTensors,
1959  subsetter,
1960  str
1961  );
1962  }
1963  }
1964 
1965 
1966  UPstream::parRun() = oldParRun;
1967 
1968 
1969  // Start sending&receiving from buffers
1970  pBufs.finishedSends();
1971 
1972 
1973  // Subset the part that stays
1974  // ~~~~~~~~~~~~~~~~~~~~~~~~~~
1975 
1976  {
1977  // Save old mesh maps before changing mesh
1978  const labelList oldFaceOwner(mesh_.faceOwner());
1979  const labelList oldFaceNeighbour(mesh_.faceNeighbour());
1980  const label oldInternalFaces = mesh_.nInternalFaces();
1981 
1982  // Remove cells.
1983  autoPtr<mapPolyMesh> subMap
1984  (
1985  doRemoveCells
1986  (
1987  select(false, distribution, Pstream::myProcNo()),
1988  oldInternalPatchi
1989  )
1990  );
1991 
1992  // Addressing from subsetted mesh
1993  subCellMap[Pstream::myProcNo()] = subMap().cellMap();
1994  subFaceMap[Pstream::myProcNo()] = renumber
1995  (
1996  repatchFaceMap,
1997  subMap().faceMap()
1998  );
1999  // Insert the sign bit from face flipping
2000  labelList& faceMap = subFaceMap[Pstream::myProcNo()];
2001  forAll(faceMap, faceI)
2002  {
2003  faceMap[faceI] += 1;
2004  }
2005  const labelHashSet& flip = subMap().flipFaceFlux();
2006  forAllConstIter(labelHashSet, flip, iter)
2007  {
2008  label faceI = iter.key();
2009  faceMap[faceI] = -faceMap[faceI];
2010  }
2011  subPointMap[Pstream::myProcNo()] = subMap().pointMap();
2012  subPatchMap[Pstream::myProcNo()] = identity(patches.size());
2013 
2014  // Initialize all addressing into current mesh
2015  constructCellMap[Pstream::myProcNo()] = identity(mesh_.nCells());
2016  constructFaceMap[Pstream::myProcNo()] = identity(mesh_.nFaces()) + 1;
2017  constructPointMap[Pstream::myProcNo()] = identity(mesh_.nPoints());
2018  constructPatchMap[Pstream::myProcNo()] = identity(patches.size());
2019 
2020  // Subset the mesh data: neighbourCell/neighbourProc
2021  // fields
2022  labelList domainSourceFace;
2023  labelList domainSourceProc;
2024  labelList domainSourcePatch;
2025  labelList domainSourceNewNbrProc;
2026 
2027  subsetBoundaryData
2028  (
2029  mesh_, // new mesh
2030  subMap().faceMap(), // from new to original mesh
2031  subMap().cellMap(),
2032 
2033  distribution, // distribution before subsetting
2034  oldFaceOwner, // owner before subsetting
2035  oldFaceNeighbour, // neighbour ,,
2036  oldInternalFaces, // nInternalFaces ,,
2037 
2038  sourceFace,
2039  sourceProc,
2040  sourcePatch,
2041  sourceNewNbrProc,
2042 
2043  domainSourceFace,
2044  domainSourceProc,
2045  domainSourcePatch,
2046  domainSourceNewNbrProc
2047  );
2048 
2049  sourceFace.transfer(domainSourceFace);
2050  sourceProc.transfer(domainSourceProc);
2051  sourcePatch.transfer(domainSourcePatch);
2052  sourceNewNbrProc.transfer(domainSourceNewNbrProc);
2053  }
2054 
2055 
2056  // Print a bit.
2057  if (debug)
2058  {
2059  Pout<< nl << "STARTING MESH:" << endl;
2060  printMeshInfo(mesh_);
2061  printFieldInfo<volScalarField>(mesh_);
2062  printFieldInfo<volVectorField>(mesh_);
2063  printFieldInfo<volSphericalTensorField>(mesh_);
2064  printFieldInfo<volSymmTensorField>(mesh_);
2065  printFieldInfo<volTensorField>(mesh_);
2066  printFieldInfo<surfaceScalarField>(mesh_);
2067  printFieldInfo<surfaceVectorField>(mesh_);
2068  printFieldInfo<surfaceSphericalTensorField>(mesh_);
2069  printFieldInfo<surfaceSymmTensorField>(mesh_);
2070  printFieldInfo<surfaceTensorField>(mesh_);
2071  Pout<< nl << endl;
2072  }
2073 
2074 
2075 
2076  // Receive and add what was sent
2077  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2078 
2079  oldParRun = UPstream::parRun();
2080  UPstream::parRun() = false;
2081 
2082  forAll(nSendCells, sendProc)
2083  {
2084  // Did processor sendProc send anything to me?
2085  if
2086  (
2087  sendProc != Pstream::myProcNo()
2088  && nSendCells[sendProc][Pstream::myProcNo()] > 0
2089  )
2090  {
2091  if (debug)
2092  {
2093  Pout<< nl
2094  << "RECEIVING FROM DOMAIN " << sendProc
2095  << " cells to receive:"
2096  << nSendCells[sendProc][Pstream::myProcNo()]
2097  << nl << endl;
2098  }
2099 
2100 
2101  // Pstream for receiving mesh and fields
2102  UIPstream str(sendProc, pBufs);
2103 
2104 
2105  // Receive from sendProc
2106  labelList domainSourceFace;
2107  labelList domainSourceProc;
2108  labelList domainSourcePatch;
2109  labelList domainSourceNewNbrProc;
2110 
2111  autoPtr<fvMesh> domainMeshPtr;
2112 
2118 
2124 
2130 
2131 
2132  // Opposite of sendMesh
2133  {
2134  domainMeshPtr = receiveMesh
2135  (
2136  sendProc,
2137  pointZoneNames,
2138  faceZoneNames,
2139  cellZoneNames,
2140 
2141  const_cast<Time&>(mesh_.time()),
2142  domainSourceFace,
2143  domainSourceProc,
2144  domainSourcePatch,
2145  domainSourceNewNbrProc,
2146  str
2147  );
2148  fvMesh& domainMesh = domainMeshPtr();
2149  // Force construction of various on mesh.
2150  //(void)domainMesh.globalData();
2151 
2152 
2153  // Receive fields. Read as single dictionary because
2154  // of problems reading consecutive fields from single stream.
2155  dictionary fieldDicts(str);
2156 
2157  // Vol fields
2158  receiveFields<volScalarField>
2159  (
2160  sendProc,
2161  volScalars,
2162  domainMesh,
2163  vsf,
2164  fieldDicts.subDict(volScalarField::typeName)
2165  );
2166  receiveFields<volVectorField>
2167  (
2168  sendProc,
2169  volVectors,
2170  domainMesh,
2171  vvf,
2172  fieldDicts.subDict(volVectorField::typeName)
2173  );
2174  receiveFields<volSphericalTensorField>
2175  (
2176  sendProc,
2177  volSphereTensors,
2178  domainMesh,
2179  vsptf,
2181  );
2182  receiveFields<volSymmTensorField>
2183  (
2184  sendProc,
2185  volSymmTensors,
2186  domainMesh,
2187  vsytf,
2189  );
2190  receiveFields<volTensorField>
2191  (
2192  sendProc,
2193  volTensors,
2194  domainMesh,
2195  vtf,
2196  fieldDicts.subDict(volTensorField::typeName)
2197  );
2198 
2199  // Surface fields
2200  receiveFields<surfaceScalarField>
2201  (
2202  sendProc,
2203  surfScalars,
2204  domainMesh,
2205  ssf,
2207  );
2208  receiveFields<surfaceVectorField>
2209  (
2210  sendProc,
2211  surfVectors,
2212  domainMesh,
2213  svf,
2215  );
2216  receiveFields<surfaceSphericalTensorField>
2217  (
2218  sendProc,
2219  surfSphereTensors,
2220  domainMesh,
2221  ssptf,
2223  );
2224  receiveFields<surfaceSymmTensorField>
2225  (
2226  sendProc,
2227  surfSymmTensors,
2228  domainMesh,
2229  ssytf,
2231  );
2232  receiveFields<surfaceTensorField>
2233  (
2234  sendProc,
2235  surfTensors,
2236  domainMesh,
2237  stf,
2239  );
2240 
2241  // Dimensioned fields
2242  receiveFields<volScalarField::Internal>
2243  (
2244  sendProc,
2245  dimScalars,
2246  domainMesh,
2247  dsf,
2248  fieldDicts.subDict
2249  (
2251  )
2252  );
2253  receiveFields<volVectorField::Internal>
2254  (
2255  sendProc,
2256  dimVectors,
2257  domainMesh,
2258  dvf,
2259  fieldDicts.subDict
2260  (
2262  )
2263  );
2264  receiveFields<volSphericalTensorField::Internal>
2265  (
2266  sendProc,
2267  dimSphereTensors,
2268  domainMesh,
2269  dstf,
2270  fieldDicts.subDict
2271  (
2273  typeName
2274  )
2275  );
2276  receiveFields<volSymmTensorField::Internal>
2277  (
2278  sendProc,
2279  dimSymmTensors,
2280  domainMesh,
2281  dsytf,
2282  fieldDicts.subDict
2283  (
2285  )
2286  );
2287  receiveFields<volTensorField::Internal>
2288  (
2289  sendProc,
2290  dimTensors,
2291  domainMesh,
2292  dtf,
2293  fieldDicts.subDict
2294  (
2296  )
2297  );
2298  }
2299  const fvMesh& domainMesh = domainMeshPtr();
2300 
2301 
2302  constructCellMap[sendProc] = identity(domainMesh.nCells());
2303  constructFaceMap[sendProc] = identity(domainMesh.nFaces()) + 1;
2304  constructPointMap[sendProc] = identity(domainMesh.nPoints());
2305  constructPatchMap[sendProc] =
2306  identity(domainMesh.boundaryMesh().size());
2307 
2308 
2309  // Print a bit.
2310  if (debug)
2311  {
2312  Pout<< nl << "RECEIVED MESH FROM:" << sendProc << endl;
2313  printMeshInfo(domainMesh);
2314  printFieldInfo<volScalarField>(domainMesh);
2315  printFieldInfo<volVectorField>(domainMesh);
2316  printFieldInfo<volSphericalTensorField>(domainMesh);
2317  printFieldInfo<volSymmTensorField>(domainMesh);
2318  printFieldInfo<volTensorField>(domainMesh);
2319  printFieldInfo<surfaceScalarField>(domainMesh);
2320  printFieldInfo<surfaceVectorField>(domainMesh);
2321  printFieldInfo<surfaceSphericalTensorField>(domainMesh);
2322  printFieldInfo<surfaceSymmTensorField>(domainMesh);
2323  printFieldInfo<surfaceTensorField>(domainMesh);
2324  }
2325 
2326 
2327  // Now this mesh we received (from sendProc) needs to be merged
2328  // with the current mesh. On the current mesh we have for all
2329  // boundaryfaces the original face and processor. See if we can
2330  // match these up to the received domainSourceFace and
2331  // domainSourceProc.
2332  labelList masterCoupledFaces;
2333  labelList slaveCoupledFaces;
2334  findCouples
2335  (
2336  mesh_,
2337 
2338  sourceFace,
2339  sourceProc,
2340  sourcePatch,
2341 
2342  sendProc,
2343  domainMesh,
2344  domainSourceFace,
2345  domainSourceProc,
2346  domainSourcePatch,
2347 
2348  masterCoupledFaces,
2349  slaveCoupledFaces
2350  );
2351 
2352  // Generate additional coupling info (points, edges) from
2353  // faces-that-match
2354  faceCoupleInfo couples
2355  (
2356  mesh_,
2357  masterCoupledFaces,
2358  domainMesh,
2359  slaveCoupledFaces,
2360  mergeTol_, // merge tolerance
2361  true, // faces align
2362  true, // couples are ordered already
2363  false
2364  );
2365 
2366 
2367  // Add domainMesh to mesh
2368  // ~~~~~~~~~~~~~~~~~~~~~~
2369 
2371  (
2372  mesh_,
2373  domainMesh,
2374  couples,
2375  false // no parallel comms
2376  );
2377 
2378  // Update mesh data: sourceFace,sourceProc for added
2379  // mesh.
2380 
2381  sourceFace = mapBoundaryData
2382  (
2383  mesh_,
2384  map(),
2385  sourceFace,
2386  domainMesh.nInternalFaces(),
2387  domainSourceFace
2388  );
2389  sourceProc = mapBoundaryData
2390  (
2391  mesh_,
2392  map(),
2393  sourceProc,
2394  domainMesh.nInternalFaces(),
2395  domainSourceProc
2396  );
2397  sourcePatch = mapBoundaryData
2398  (
2399  mesh_,
2400  map(),
2401  sourcePatch,
2402  domainMesh.nInternalFaces(),
2403  domainSourcePatch
2404  );
2405  sourceNewNbrProc = mapBoundaryData
2406  (
2407  mesh_,
2408  map(),
2409  sourceNewNbrProc,
2410  domainMesh.nInternalFaces(),
2411  domainSourceNewNbrProc
2412  );
2413 
2414  // Update all addressing so xxProcAddressing points to correct
2415  // item in masterMesh.
2416  const labelList& oldCellMap = map().oldCellMap();
2417  const labelList& oldFaceMap = map().oldFaceMap();
2418  const labelList& oldPointMap = map().oldPointMap();
2419  const labelList& oldPatchMap = map().oldPatchMap();
2420 
2421  //Note: old mesh faces never flipped!
2422  forAll(constructPatchMap, proci)
2423  {
2424  if (proci != sendProc && constructPatchMap[proci].size())
2425  {
2426  // Processor already in mesh (either myProcNo or received)
2427  inplaceRenumber(oldCellMap, constructCellMap[proci]);
2428  inplaceRenumberWithFlip
2429  (
2430  oldFaceMap,
2431  false,
2432  true,
2433  constructFaceMap[proci]
2434  );
2435  inplaceRenumber(oldPointMap, constructPointMap[proci]);
2436  inplaceRenumber(oldPatchMap, constructPatchMap[proci]);
2437  }
2438  }
2439 
2440 
2441  labelHashSet flippedAddedFaces;
2442  {
2443  // Find out if any faces of domain mesh were flipped (boundary
2444  // faces becoming internal)
2445  label nBnd = domainMesh.nFaces()-domainMesh.nInternalFaces();
2446  flippedAddedFaces.resize(nBnd/4);
2447 
2448  for
2449  (
2450  label domainFaceI = domainMesh.nInternalFaces();
2451  domainFaceI < domainMesh.nFaces();
2452  domainFaceI++
2453  )
2454  {
2455  label newFaceI = map().addedFaceMap()[domainFaceI];
2456  label newCellI = mesh_.faceOwner()[newFaceI];
2457 
2458  label domainCellI = domainMesh.faceOwner()[domainFaceI];
2459 
2460  if (newCellI != map().addedCellMap()[domainCellI])
2461  {
2462  flippedAddedFaces.insert(domainFaceI);
2463  }
2464  }
2465  }
2466 
2467 
2468  // Added processor
2469  inplaceRenumber(map().addedCellMap(), constructCellMap[sendProc]);
2470  // Add flip
2471  forAllConstIter(labelHashSet, flippedAddedFaces, iter)
2472  {
2473  label domainFaceI = iter.key();
2474  label& val = constructFaceMap[sendProc][domainFaceI];
2475  val = -val;
2476  }
2477  inplaceRenumberWithFlip
2478  (
2479  map().addedFaceMap(),
2480  false,
2481  true, // constructFaceMap has flip sign
2482  constructFaceMap[sendProc]
2483  );
2484  inplaceRenumber(map().addedPointMap(), constructPointMap[sendProc]);
2485  inplaceRenumber(map().addedPatchMap(), constructPatchMap[sendProc]);
2486 
2487  if (debug)
2488  {
2489  Pout<< nl << "MERGED MESH FROM:" << sendProc << endl;
2490  printMeshInfo(mesh_);
2491  printFieldInfo<volScalarField>(mesh_);
2492  printFieldInfo<volVectorField>(mesh_);
2493  printFieldInfo<volSphericalTensorField>(mesh_);
2494  printFieldInfo<volSymmTensorField>(mesh_);
2495  printFieldInfo<volTensorField>(mesh_);
2496  printFieldInfo<surfaceScalarField>(mesh_);
2497  printFieldInfo<surfaceVectorField>(mesh_);
2498  printFieldInfo<surfaceSphericalTensorField>(mesh_);
2499  printFieldInfo<surfaceSymmTensorField>(mesh_);
2500  printFieldInfo<surfaceTensorField>(mesh_);
2501  Pout<< nl << endl;
2502  }
2503  }
2504  }
2505 
2506  UPstream::parRun() = oldParRun;
2507 
2508  // Print a bit.
2509  if (debug)
2510  {
2511  Pout<< nl << "REDISTRIBUTED MESH:" << endl;
2512  printMeshInfo(mesh_);
2513  printFieldInfo<volScalarField>(mesh_);
2514  printFieldInfo<volVectorField>(mesh_);
2515  printFieldInfo<volSphericalTensorField>(mesh_);
2516  printFieldInfo<volSymmTensorField>(mesh_);
2517  printFieldInfo<volTensorField>(mesh_);
2518  printFieldInfo<surfaceScalarField>(mesh_);
2519  printFieldInfo<surfaceVectorField>(mesh_);
2520  printFieldInfo<surfaceSphericalTensorField>(mesh_);
2521  printFieldInfo<surfaceSymmTensorField>(mesh_);
2522  printFieldInfo<surfaceTensorField>(mesh_);
2523  Pout<< nl << endl;
2524  }
2525 
2526 
2527 
2528  // Add processorPatches
2529  // ~~~~~~~~~~~~~~~~~~~~
2530 
2531  // Per neighbour processor, per originating patch, the patchID
2532  // For faces resulting from internal faces or normal processor patches
2533  // the originating patch is -1. For cyclics this is the cyclic patchID.
2534  List<Map<label>> procPatchID;
2535 
2536  // Add processor and processorCyclic patches.
2537  addProcPatches(sourceNewNbrProc, sourcePatch, procPatchID);
2538 
2539  // Put faces into correct patch. Note that we now have proper
2540  // processorPolyPatches again so repatching will take care of coupled face
2541  // ordering.
2542 
2543  // Get boundary faces to be repatched. Is -1 or new patchID
2544  labelList newPatchID
2545  (
2546  getBoundaryPatch
2547  (
2548  sourceNewNbrProc,
2549  sourcePatch,
2550  procPatchID
2551  )
2552  );
2553 
2554  // Change patches. Since this might change ordering of coupled faces
2555  // we also need to adapt our constructMaps.
2556  // NOTE: there is one very particular problem with this structure.
2557  // We first create the processor patches and use these to merge out
2558  // shared points (see mergeSharedPoints below). So temporarily points
2559  // and edges do not match!
2560  repatch(newPatchID, constructFaceMap);
2561 
2562  // See if any geometrically shared points need to be merged. Note: does
2563  // parallel comms. After this points and edges should again be consistent.
2564  mergeSharedPoints(constructPointMap);
2565 
2566  // Bit of hack: processorFvPatchField does not get reset since created
2567  // from nothing so explicitly reset.
2568  initPatchFields<volScalarField, processorFvPatchField<scalar>>
2569  (
2570  Zero
2571  );
2572  initPatchFields<volVectorField, processorFvPatchField<vector>>
2573  (
2574  Zero
2575  );
2576  initPatchFields
2577  <
2580  >
2581  (
2582  Zero
2583  );
2584  initPatchFields<volSymmTensorField, processorFvPatchField<symmTensor>>
2585  (
2586  Zero
2587  );
2588  initPatchFields<volTensorField, processorFvPatchField<tensor>>
2589  (
2590  Zero
2591  );
2592 
2593 
2594  mesh_.setInstance(mesh_.time().timeName());
2595 
2596 
2597  // Print a bit
2598  if (debug)
2599  {
2600  Pout<< nl << "FINAL MESH:" << endl;
2601  printMeshInfo(mesh_);
2602  printFieldInfo<volScalarField>(mesh_);
2603  printFieldInfo<volVectorField>(mesh_);
2604  printFieldInfo<volSphericalTensorField>(mesh_);
2605  printFieldInfo<volSymmTensorField>(mesh_);
2606  printFieldInfo<volTensorField>(mesh_);
2607  printFieldInfo<surfaceScalarField>(mesh_);
2608  printFieldInfo<surfaceVectorField>(mesh_);
2609  printFieldInfo<surfaceSphericalTensorField>(mesh_);
2610  printFieldInfo<surfaceSymmTensorField>(mesh_);
2611  printFieldInfo<surfaceTensorField>(mesh_);
2612  Pout<< nl << endl;
2613  }
2614 
2615  // Collect all maps and return
2617  (
2619  (
2620  mesh_,
2621 
2622  nOldPoints,
2623  nOldFaces,
2624  nOldCells,
2625  oldPatchStarts.xfer(),
2626  oldPatchNMeshPoints.xfer(),
2627 
2628  subPointMap.xfer(),
2629  subFaceMap.xfer(),
2630  subCellMap.xfer(),
2631  subPatchMap.xfer(),
2632 
2633  constructPointMap.xfer(),
2634  constructFaceMap.xfer(),
2635  constructCellMap.xfer(),
2636  constructPatchMap.xfer(),
2637 
2638  true, // subFaceMap has flip
2639  true // constructFaceMap has flip
2640  )
2641  );
2642 }
2643 
2644 
2645 // ************************************************************************* //
Foam::surfaceFields.
ListType renumber(const labelUList &oldToNew, const ListType &)
Renumber the values (not the indices) of a list.
dimensionedScalar sign(const dimensionedScalar &ds)
const labelList & faceFlipMap() const
Return face map with sign to encode flipped faces.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
void clearOut()
Clear all geometry and addressing.
Definition: fvMesh.C:230
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
A HashTable with keys but without contents.
Definition: HashSet.H:59
dictionary dict
static void reorderPatches(fvMesh &, const labelList &oldToNew, const label nPatches, const bool validBoundary)
Reorder and remove trailing patches. If validBoundary call is parallel.
Definition: fvMeshTools.C:324
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
void setRefinement(const labelList &cellsToRemove, const labelList &facesToExpose, const labelList &patchIDs, polyTopoChange &) const
Play commands into polyTopoChange to remove cells.
Definition: removeCells.C:179
void finishedSends(const bool block=true)
Mark all sends as having been done. This will start receives.
const labelList & pointMap() const
Return point map.
bool checkParallelSync(const bool report=false) const
Check whether all procs have all patches and in same order. Return.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
wordList names() const
Return a list of patch names.
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:253
Class describing modification of a face.
void resetMotion() const
Reset motion.
Definition: polyMesh.C:1133
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Xfer< List< T > > xfer()
Transfer contents to the Xfer container.
Definition: ListI.H:90
static const char *const typeName
Definition: Field.H:94
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
int neighbProcNo() const
Return neigbour processor number.
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:633
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:417
Given list of cells to remove insert all the topology changes.
Definition: removeCells.H:59
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:306
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: UList.H:440
Class containing mesh-to-mesh mapping information after a mesh addition where we add a mesh (&#39;added m...
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
Less function class that can be used for sorting processor patches.
Class containing mesh-to-mesh mapping information after a mesh distribution where we send parts of me...
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:715
patches[0]
const labelList & oldFaceMap() const
const vectorField & faceCentres() const
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
GeometricField< sphericalTensor, fvPatchField, volMesh > volSphericalTensorField
Definition: volFieldsFwd.H:57
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
void addFvPatches(const List< polyPatch * > &, const bool validBoundary=true)
Add boundary patches. Constructor helper.
Definition: fvMesh.C:460
Combination-Reduction operation for a parallel run. The information from all nodes is collected on th...
virtual transformType transform() const
Type of transform.
const fvMesh & subMesh() const
Return reference to subset mesh.
const labelList & addedFaceMap() const
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
label nCells() const
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
Various functions to operate on Lists.
Container for information needed to couple to meshes. When constructed from two meshes and a geometri...
Neighbour processor patch.
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:305
Xfer< HashTable< T, Key, Hash > > xfer()
Transfer contents to the Xfer container.
Definition: HashTableI.H:102
Input inter-processor communications stream operating on external buffer.
Definition: UIPstream.H:53
const word & name() const
Return name.
Definition: zone.H:150
A List obtained as a section of another List.
Definition: SubList.H:53
label findZoneID(const word &zoneName) const
Find zone index given a name.
Definition: ZoneMesh.C:341
const labelList & faceMap() const
Return face map.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
bool operator()(const label a, const label b)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
static void mergePoints(const polyMesh &, const Map< label > &pointToMaster, polyTopoChange &meshMod)
Helper: Merge points.
Xfer< T > xferMove(T &)
Construct by transferring the contents of the arg.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:310
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:32
Cyclic plane patch.
wordList names() const
Return a list of zone names.
Definition: ZoneMesh.C:256
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:91
const labelList & patchMap() const
Return patch map.
labelList getExposedFaces(const labelList &cellsToRemove) const
Get labels of exposed faces.
Definition: removeCells.C:73
wordList names() const
Return the list of names of the IOobjects.
const word & name() const
Return name.
An STL-conforming hash table.
Definition: HashTable.H:61
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
List< Key > toc() const
Return the table of contents.
static autoPtr< mapAddedPolyMesh > add(fvMesh &mesh0, const fvMesh &mesh1, const faceCoupleInfo &coupleInfo, const bool validBoundary=true)
Inplace add mesh to fvMesh. Maps all stored fields. Returns map.
Definition: fvMeshAdder.C:71
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
Foam::polyBoundaryMesh.
A packed storage unstructured matrix of objects of type <T> using an offset table for access...
void setSize(const label mRows)
Reset size of CompactListList.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
Output inter-processor communications stream operating on external buffer.
Definition: UOPstream.H:54
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
static const char nl
Definition: Ostream.H:262
defineTypeNameAndDebug(combustionModel, 0)
autoPtr< mapDistributePolyMesh > distribute(const labelList &dist)
Send cells to neighbours according to distribution.
Post-processing mesh subset tool. Given the original mesh and the list of selected cells...
Definition: fvMeshSubset.H:73
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
Buffers for inter-processor communications streams (UOPstream, UIPstream).
const pointZoneMesh & pointZones() const
Return point zone mesh.
Definition: polyMesh.H:457
A subset of mesh cells.
Definition: cellZone.H:61
label nOldInternalFaces() const
Number of old internal faces.
void setInstance(const fileName &)
Set the instance for mesh files.
Definition: polyMeshIO.C:32
This boundary condition enables processor communication across cyclic patches.
const labelList & cellMap() const
Return cell map.
label nFaces() const
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
void setSize(const label)
Reset size of List.
Definition: List.C:295
static void printMeshInfo(const fvMesh &)
Print some info on mesh.
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:393
const cellZoneMesh & cellZones() const
Return cell zone mesh.
Definition: polyMesh.H:469
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:399
label patchi
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1023
A subset of mesh points. The labels of points in the zone can be obtained from the addressing() list...
Definition: pointZone.H:62
virtual bool owner() const
Does the processor own the patch ?
const labelUList & faceCells() const
Return face-cell addressing.
Definition: polyPatch.C:315
label newPointi
Definition: readKivaGrid.H:501
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:62
Foam::fvBoundaryMesh.
static void printCoupleInfo(const primitiveMesh &, const labelList &, const labelList &, const labelList &, const labelList &)
Print some info on coupling data.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
A List with indirect addressing.
Definition: fvMatrix.H:106
Direct mesh changes based on v1.3 polyTopoChange syntax.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
lessProcPatches(const labelList &nbrProc, const labelList &referPatchID)
label nPoints() const
virtual bool owner() const
Does this side own the patch ?
label n
void resize(const label newSize)
Resize the hash table for efficiency.
Definition: HashTable.C:428
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:463
void addZones(const List< pointZone * > &pz, const List< faceZone * > &fz, const List< cellZone * > &cz)
Add mesh zones.
Definition: polyMesh.C:922
static autoPtr< polyPatch > New(const word &patchType, const word &name, const label size, const label start, const label index, const polyBoundaryMesh &bm)
Return a pointer to a new patch created on freestore from.
Definition: polyPatchNew.C:32
This boundary condition enables processor communication across patches.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:53
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:430
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1017
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
static labelList countCells(const labelList &)
Helper function: count cells per processor in wanted distribution.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
const polyBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
Definition: polyPatch.C:278
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:545
label nInternalFaces() const
bool found
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:365
void setLargeCellSubset(const labelList &region, const label currentRegion, const label patchID=-1, const bool syncCouples=true)
Set the subset from all cells with region == currentRegion.
Definition: fvMeshSubset.C:795
static Map< label > findSharedPoints(const polyMesh &, const scalar mergeTol)
Find topologically and geometrically shared points.
const word & name() const
Return name.
Definition: IOobject.H:260
Namespace for OpenFOAM.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.