FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 
27 #include "fvMeshDistribute.H"
37 #include <OpenFOAM/mergePoints.H>
40 #include <OpenFOAM/syncTools.H>
41 
42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43 
45 
46 
47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 
49 Foam::labelList Foam::fvMeshDistribute::select
50 (
51  const bool selectEqual,
52  const labelList& values,
53  const label value
54 )
55 {
56  label n = 0;
57 
58  forAll(values, i)
59  {
60  if (selectEqual == (values[i] == value))
61  {
62  n++;
63  }
64  }
65 
66  labelList indices(n);
67  n = 0;
68 
69  forAll(values, i)
70  {
71  if (selectEqual == (values[i] == value))
72  {
73  indices[n++] = i;
74  }
75  }
76  return indices;
77 }
78 
79 
80 // Check all procs have same names and in exactly same order.
81 void Foam::fvMeshDistribute::checkEqualWordList
82 (
83  const string& msg,
84  const wordList& lst
85 )
86 {
87  List<wordList> allNames(Pstream::nProcs());
88  allNames[Pstream::myProcNo()] = lst;
89  Pstream::gatherList(allNames);
90  Pstream::scatterList(allNames);
91 
92  for (label procI = 1; procI < Pstream::nProcs(); procI++)
93  {
94  if (allNames[procI] != allNames[0])
95  {
96  FatalErrorIn("fvMeshDistribute::checkEqualWordList(..)")
97  << "When checking for equal " << msg.c_str() << " :" << endl
98  << "processor0 has:" << allNames[0] << endl
99  << "processor" << procI << " has:" << allNames[procI] << endl
100  << msg.c_str() << " need to be synchronised on all processors."
101  << exit(FatalError);
102  }
103  }
104 }
105 
106 
107 Foam::wordList Foam::fvMeshDistribute::mergeWordList(const wordList& procNames)
108 {
109  List<wordList> allNames(Pstream::nProcs());
110  allNames[Pstream::myProcNo()] = procNames;
111  Pstream::gatherList(allNames);
112  Pstream::scatterList(allNames);
113 
114  HashSet<word> mergedNames;
115  forAll(allNames, procI)
116  {
117  forAll(allNames[procI], i)
118  {
119  mergedNames.insert(allNames[procI][i]);
120  }
121  }
122  return mergedNames.toc();
123 }
124 
125 
126 // Print some info on mesh.
128 {
129  Pout<< "Primitives:" << nl
130  << " points :" << mesh.nPoints() << nl
131  << " internalFaces:" << mesh.nInternalFaces() << nl
132  << " faces :" << mesh.nFaces() << nl
133  << " cells :" << mesh.nCells() << nl;
134 
135  const fvBoundaryMesh& patches = mesh.boundary();
136 
137  Pout<< "Patches:" << endl;
138  forAll(patches, patchI)
139  {
140  const polyPatch& pp = patches[patchI].patch();
141 
142  Pout<< " " << patchI << " name:" << pp.name()
143  << " size:" << pp.size()
144  << " start:" << pp.start()
145  << " type:" << pp.type()
146  << endl;
147  }
148 
149  if (mesh.pointZones().size())
150  {
151  Pout<< "PointZones:" << endl;
152  forAll(mesh.pointZones(), zoneI)
153  {
154  const pointZone& pz = mesh.pointZones()[zoneI];
155  Pout<< " " << zoneI << " name:" << pz.name()
156  << " size:" << pz.size()
157  << endl;
158  }
159  }
160  if (mesh.faceZones().size())
161  {
162  Pout<< "FaceZones:" << endl;
163  forAll(mesh.faceZones(), zoneI)
164  {
165  const faceZone& fz = mesh.faceZones()[zoneI];
166  Pout<< " " << zoneI << " name:" << fz.name()
167  << " size:" << fz.size()
168  << endl;
169  }
170  }
171  if (mesh.cellZones().size())
172  {
173  Pout<< "CellZones:" << endl;
174  forAll(mesh.cellZones(), zoneI)
175  {
176  const cellZone& cz = mesh.cellZones()[zoneI];
177  Pout<< " " << zoneI << " name:" << cz.name()
178  << " size:" << cz.size()
179  << endl;
180  }
181  }
182 }
183 
184 
186 (
187  const primitiveMesh& mesh,
188  const labelList& sourceFace,
189  const labelList& sourceProc,
190  const labelList& sourceNewProc
191 )
192 {
193  Pout<< nl
194  << "Current coupling info:"
195  << endl;
196 
197  forAll(sourceFace, bFaceI)
198  {
199  label meshFaceI = mesh.nInternalFaces() + bFaceI;
200 
201  Pout<< " meshFace:" << meshFaceI
202  << " fc:" << mesh.faceCentres()[meshFaceI]
203  << " connects to proc:" << sourceProc[bFaceI]
204  << "/face:" << sourceFace[bFaceI]
205  << " which will move to proc:" << sourceNewProc[bFaceI]
206  << endl;
207  }
208 }
209 
210 
211 // Finds (non-empty) patch that exposed internal and proc faces can be put into.
212 Foam::label Foam::fvMeshDistribute::findNonEmptyPatch() const
213 {
214  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
215 
216  label nonEmptyPatchI = -1;
217 
218  forAllReverse(patches, patchI)
219  {
220  const polyPatch& pp = patches[patchI];
221 
222  if (!isA<emptyPolyPatch>(pp) && !pp.coupled())
223  {
224  nonEmptyPatchI = patchI;
225  break;
226  }
227  }
228 
229  if (nonEmptyPatchI == -1)
230  {
231  FatalErrorIn("fvMeshDistribute::findNonEmptyPatch() const")
232  << "Cannot find a patch which is neither of type empty nor"
233  << " coupled in patches " << patches.names() << endl
234  << "There has to be at least one such patch for"
235  << " distribution to work" << abort(FatalError);
236  }
237 
238  if (debug)
239  {
240  Pout<< "findNonEmptyPatch : using patch " << nonEmptyPatchI
241  << " name:" << patches[nonEmptyPatchI].name()
242  << " type:" << patches[nonEmptyPatchI].type()
243  << " to put exposed faces into." << endl;
244  }
245 
246 
247  // Do additional test for processor patches intermingled with non-proc
248  // patches.
249  label procPatchI = -1;
250 
251  forAll(patches, patchI)
252  {
253  if (isA<processorPolyPatch>(patches[patchI]))
254  {
255  procPatchI = patchI;
256  }
257  else if (procPatchI != -1)
258  {
259  FatalErrorIn("fvMeshDistribute::findNonEmptyPatch() const")
260  << "Processor patches should be at end of patch list."
261  << endl
262  << "Have processor patch " << procPatchI
263  << " followed by non-processor patch " << patchI
264  << " in patches " << patches.names()
265  << abort(FatalError);
266  }
267  }
268 
269  return nonEmptyPatchI;
270 }
271 
272 
273 // Appends processorPolyPatch. Returns patchID.
274 Foam::label Foam::fvMeshDistribute::addProcPatch
275 (
276  const word& patchName,
277  const label nbrProc
278 )
279 {
280  // Clear local fields and e.g. polyMesh globalMeshData.
281  mesh_.clearOut();
282 
283 
284  polyBoundaryMesh& polyPatches =
285  const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh());
286  fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh_.boundary());
287 
288  if (polyPatches.findPatchID(patchName) != -1)
289  {
290  FatalErrorIn("fvMeshDistribute::addProcPatch(const word&, const label)")
291  << "Cannot create patch " << patchName << " since already exists."
292  << nl
293  << "Current patch names:" << polyPatches.names()
294  << exit(FatalError);
295  }
296 
297 
298 
299  // Add the patch
300  // ~~~~~~~~~~~~~
301 
302  label sz = polyPatches.size();
303 
304  // Add polyPatch
305  polyPatches.setSize(sz+1);
306  polyPatches.set
307  (
308  sz,
309  new processorPolyPatch
310  (
311  patchName,
312  0, // size
313  mesh_.nFaces(),
314  sz,
315  mesh_.boundaryMesh(),
317  nbrProc
318  )
319  );
320  fvPatches.setSize(sz+1);
321  fvPatches.set
322  (
323  sz,
325  (
326  polyPatches[sz], // point to newly added polyPatch
327  mesh_.boundary()
328  )
329  );
330 
331  return sz;
332 }
333 
334 
335 // Deletes last patch
336 void Foam::fvMeshDistribute::deleteTrailingPatch()
337 {
338  // Clear local fields and e.g. polyMesh globalMeshData.
339  mesh_.clearOut();
340 
341  polyBoundaryMesh& polyPatches =
342  const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh());
343  fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh_.boundary());
344 
345  if (polyPatches.empty())
346  {
347  FatalErrorIn("fvMeshDistribute::deleteTrailingPatch(fvMesh&)")
348  << "No patches in mesh"
349  << abort(FatalError);
350  }
351 
352  label sz = polyPatches.size();
353 
354  label nFaces = polyPatches[sz-1].size();
355 
356  // Remove last polyPatch
357  if (debug)
358  {
359  Pout<< "deleteTrailingPatch : Removing patch " << sz-1
360  << " : " << polyPatches[sz-1].name() << " size:" << nFaces << endl;
361  }
362 
363  if (nFaces)
364  {
365  FatalErrorIn("fvMeshDistribute::deleteTrailingPatch()")
366  << "There are still " << nFaces << " faces in patch to be deleted "
367  << sz-1 << ' ' << polyPatches[sz-1].name()
368  << abort(FatalError);
369  }
370 
371  // Remove actual patch
372  polyPatches.setSize(sz-1);
373  fvPatches.setSize(sz-1);
374 }
375 
376 
377 // Delete all processor patches. Move any processor faces into the last
378 // non-processor patch.
379 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::deleteProcPatches
380 (
381  const label destinationPatch
382 )
383 {
384  // New patchID per boundary faces to be repatched. Is -1 (no change)
385  // or new patchID
386  labelList newPatchID(mesh_.nFaces() - mesh_.nInternalFaces(), -1);
387 
388  label nProcPatches = 0;
389 
390  forAll(mesh_.boundaryMesh(), patchI)
391  {
392  const polyPatch& pp = mesh_.boundaryMesh()[patchI];
393 
394  if (isA<processorPolyPatch>(pp))
395  {
396  if (debug)
397  {
398  Pout<< "Moving all faces of patch " << pp.name()
399  << " into patch " << destinationPatch
400  << endl;
401  }
402 
403  label offset = pp.start() - mesh_.nInternalFaces();
404 
405  forAll(pp, i)
406  {
407  newPatchID[offset+i] = destinationPatch;
408  }
409 
410  nProcPatches++;
411  }
412  }
413 
414  // Note: order of boundary faces been kept the same since the
415  // destinationPatch is at the end and we have visited the patches in
416  // incremental order.
417  labelListList dummyFaceMaps;
418  autoPtr<mapPolyMesh> map = repatch(newPatchID, dummyFaceMaps);
419 
420 
421  // Delete (now empty) processor patches.
422  forAllReverse(mesh_.boundaryMesh(), patchI)
423  {
424  const polyPatch& pp = mesh_.boundaryMesh()[patchI];
425 
426  if (isA<processorPolyPatch>(pp))
427  {
428  deleteTrailingPatch();
429  deleteTrailingPatchFields<volScalarField>();
430  deleteTrailingPatchFields<volVectorField>();
431  deleteTrailingPatchFields<volSphericalTensorField>();
432  deleteTrailingPatchFields<volSymmTensorField>();
433  deleteTrailingPatchFields<volTensorField>();
434 
435  deleteTrailingPatchFields<surfaceScalarField>();
436  deleteTrailingPatchFields<surfaceVectorField>();
437  deleteTrailingPatchFields<surfaceSphericalTensorField>();
438  deleteTrailingPatchFields<surfaceSymmTensorField>();
439  deleteTrailingPatchFields<surfaceTensorField>();
440  }
441  }
442 
443  return map;
444 }
445 
446 
447 // Repatch the mesh.
448 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::repatch
449 (
450  const labelList& newPatchID, // per boundary face -1 or new patchID
451  labelListList& constructFaceMap
452 )
453 {
454  polyTopoChange meshMod(mesh_);
455 
456  forAll(newPatchID, bFaceI)
457  {
458  if (newPatchID[bFaceI] != -1)
459  {
460  label faceI = mesh_.nInternalFaces() + bFaceI;
461 
462  label zoneID = mesh_.faceZones().whichZone(faceI);
463  bool zoneFlip = false;
464 
465  if (zoneID >= 0)
466  {
467  const faceZone& fZone = mesh_.faceZones()[zoneID];
468  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
469  }
470 
471  meshMod.setAction
472  (
473  polyModifyFace
474  (
475  mesh_.faces()[faceI], // modified face
476  faceI, // label of face
477  mesh_.faceOwner()[faceI], // owner
478  -1, // neighbour
479  false, // face flip
480  newPatchID[bFaceI], // patch for face
481  false, // remove from zone
482  zoneID, // zone for face
483  zoneFlip // face flip in zone
484  )
485  );
486  }
487  }
488 
489 
490  // Do mapping of fields from one patchField to the other ourselves since
491  // is currently not supported by updateMesh.
492 
493  // Store boundary fields (we only do this for surfaceFields)
494  PtrList<FieldField<fvsPatchField, scalar> > sFlds;
495  saveBoundaryFields<scalar, surfaceMesh>(sFlds);
496  PtrList<FieldField<fvsPatchField, vector> > vFlds;
497  saveBoundaryFields<vector, surfaceMesh>(vFlds);
498  PtrList<FieldField<fvsPatchField, sphericalTensor> > sptFlds;
499  saveBoundaryFields<sphericalTensor, surfaceMesh>(sptFlds);
500  PtrList<FieldField<fvsPatchField, symmTensor> > sytFlds;
501  saveBoundaryFields<symmTensor, surfaceMesh>(sytFlds);
502  PtrList<FieldField<fvsPatchField, tensor> > tFlds;
503  saveBoundaryFields<tensor, surfaceMesh>(tFlds);
504 
505  // Change the mesh (no inflation). Note: parallel comms allowed.
506  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
507 
508  // Update fields. No inflation, parallel sync.
509  mesh_.updateMesh(map);
510 
511  // Map patch fields using stored boundary fields. Note: assumes order
512  // of fields has not changed in object registry!
513  mapBoundaryFields<scalar, surfaceMesh>(map, sFlds);
514  mapBoundaryFields<vector, surfaceMesh>(map, vFlds);
515  mapBoundaryFields<sphericalTensor, surfaceMesh>(map, sptFlds);
516  mapBoundaryFields<symmTensor, surfaceMesh>(map, sytFlds);
517  mapBoundaryFields<tensor, surfaceMesh>(map, tFlds);
518 
519 
520  // Move mesh (since morphing does not do this)
521  if (map().hasMotionPoints())
522  {
523  mesh_.movePoints(map().preMotionPoints());
524  }
525 
526  // Adapt constructMaps.
527 
528  if (debug)
529  {
530  label index = findIndex(map().reverseFaceMap(), -1);
531 
532  if (index != -1)
533  {
535  (
536  "fvMeshDistribute::repatch(const labelList&, labelListList&)"
537  ) << "reverseFaceMap contains -1 at index:"
538  << index << endl
539  << "This means that the repatch operation was not just"
540  << " a shuffle?" << abort(FatalError);
541  }
542  }
543 
544  forAll(constructFaceMap, procI)
545  {
546  inplaceRenumber(map().reverseFaceMap(), constructFaceMap[procI]);
547  }
548 
549 
550  return map;
551 }
552 
553 
554 // Detect shared points. Need processor patches to be present.
555 // Background: when adding bits of mesh one can get points which
556 // share the same position but are only detectable to be topologically
557 // the same point when doing parallel analysis. This routine will
558 // merge those points.
559 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
560 (
561  labelListList& constructPointMap
562 )
563 {
564  // Find out which sets of points get merged and create a map from
565  // mesh point to unique point.
566  Map<label> pointToMaster
567  (
569  (
570  mesh_,
571  mergeTol_
572  )
573  );
574 
575  if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
576  {
577  return autoPtr<mapPolyMesh>(NULL);
578  }
579 
580  polyTopoChange meshMod(mesh_);
581 
582  fvMeshAdder::mergePoints(mesh_, pointToMaster, meshMod);
583 
584  // Change the mesh (no inflation). Note: parallel comms allowed.
585  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
586 
587  // Update fields. No inflation, parallel sync.
588  mesh_.updateMesh(map);
589 
590  // Adapt constructMaps for merged points.
591  forAll(constructPointMap, procI)
592  {
593  labelList& constructMap = constructPointMap[procI];
594 
595  forAll(constructMap, i)
596  {
597  label oldPointI = constructMap[i];
598 
599  label newPointI = map().reversePointMap()[oldPointI];
600 
601  if (newPointI < -1)
602  {
603  constructMap[i] = -newPointI-2;
604  }
605  else if (newPointI >= 0)
606  {
607  constructMap[i] = newPointI;
608  }
609  else
610  {
611  FatalErrorIn("fvMeshDistribute::mergeSharedPoints()")
612  << "Problem. oldPointI:" << oldPointI
613  << " newPointI:" << newPointI << abort(FatalError);
614  }
615  }
616  }
617  return map;
618 }
619 
620 
621 // Construct the local environment of all boundary faces.
622 void Foam::fvMeshDistribute::getNeighbourData
623 (
624  const labelList& distribution,
625  labelList& sourceFace,
626  labelList& sourceProc,
627  labelList& sourceNewProc
628 ) const
629 {
630  label nBnd = mesh_.nFaces() - mesh_.nInternalFaces();
631  sourceFace.setSize(nBnd);
632  sourceProc.setSize(nBnd);
633  sourceNewProc.setSize(nBnd);
634 
635  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
636 
637  // Get neighbouring meshFace labels and new processor of coupled boundaries.
638  labelList nbrFaces(nBnd, -1);
639  labelList nbrNewProc(nBnd, -1);
640 
641  forAll(patches, patchI)
642  {
643  const polyPatch& pp = patches[patchI];
644 
645  if (isA<processorPolyPatch>(pp))
646  {
647  label offset = pp.start() - mesh_.nInternalFaces();
648 
649  // Mesh labels of faces on this side
650  forAll(pp, i)
651  {
652  label bndI = offset + i;
653  nbrFaces[bndI] = pp.start()+i;
654  }
655 
656  // Which processor they will end up on
657  SubList<label>(nbrNewProc, pp.size(), offset).assign
658  (
659  UIndirectList<label>(distribution, pp.faceCells())()
660  );
661  }
662  }
663 
664 
665  // Exchange the boundary data
666  syncTools::swapBoundaryFaceList(mesh_, nbrFaces, false);
667  syncTools::swapBoundaryFaceList(mesh_, nbrNewProc, false);
668 
669 
670  forAll(patches, patchI)
671  {
672  const polyPatch& pp = patches[patchI];
673  label offset = pp.start() - mesh_.nInternalFaces();
674 
675  if (isA<processorPolyPatch>(pp))
676  {
677  const processorPolyPatch& procPatch =
678  refCast<const processorPolyPatch>(pp);
679 
680  // Check which of the two faces we store.
681 
682  if (Pstream::myProcNo() < procPatch.neighbProcNo())
683  {
684  // Use my local face labels
685  forAll(pp, i)
686  {
687  label bndI = offset + i;
688  sourceFace[bndI] = pp.start()+i;
689  sourceProc[bndI] = Pstream::myProcNo();
690  sourceNewProc[bndI] = nbrNewProc[bndI];
691  }
692  }
693  else
694  {
695  // Use my neighbours face labels
696  forAll(pp, i)
697  {
698  label bndI = offset + i;
699  sourceFace[bndI] = nbrFaces[bndI];
700  sourceProc[bndI] = procPatch.neighbProcNo();
701  sourceNewProc[bndI] = nbrNewProc[bndI];
702  }
703  }
704  }
705  else
706  {
707  // Normal (physical) boundary
708  forAll(pp, i)
709  {
710  label bndI = offset + i;
711  sourceFace[bndI] = patchI;
712  sourceProc[bndI] = -1;
713  sourceNewProc[bndI] = -1;
714  }
715  }
716  }
717 }
718 
719 
720 // Subset the neighbourCell/neighbourProc fields
721 void Foam::fvMeshDistribute::subsetBoundaryData
722 (
723  const fvMesh& mesh,
724  const labelList& faceMap,
725  const labelList& cellMap,
726 
727  const labelList& oldDistribution,
728  const labelList& oldFaceOwner,
729  const labelList& oldFaceNeighbour,
730  const label oldInternalFaces,
731 
732  const labelList& sourceFace,
733  const labelList& sourceProc,
734  const labelList& sourceNewProc,
735 
736  labelList& subFace,
737  labelList& subProc,
738  labelList& subNewProc
739 )
740 {
741  subFace.setSize(mesh.nFaces() - mesh.nInternalFaces());
742  subProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
743  subNewProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
744 
745  forAll(subFace, newBFaceI)
746  {
747  label newFaceI = newBFaceI + mesh.nInternalFaces();
748 
749  label oldFaceI = faceMap[newFaceI];
750 
751  // Was oldFaceI internal face? If so which side did we get.
752  if (oldFaceI < oldInternalFaces)
753  {
754  subFace[newBFaceI] = oldFaceI;
755  subProc[newBFaceI] = Pstream::myProcNo();
756 
757  label oldOwn = oldFaceOwner[oldFaceI];
758  label oldNei = oldFaceNeighbour[oldFaceI];
759 
760  if (oldOwn == cellMap[mesh.faceOwner()[newFaceI]])
761  {
762  // We kept the owner side. Where does the neighbour move to?
763  subNewProc[newBFaceI] = oldDistribution[oldNei];
764  }
765  else
766  {
767  // We kept the neighbour side.
768  subNewProc[newBFaceI] = oldDistribution[oldOwn];
769  }
770  }
771  else
772  {
773  // Was boundary face. Take over boundary information
774  label oldBFaceI = oldFaceI - oldInternalFaces;
775 
776  subFace[newBFaceI] = sourceFace[oldBFaceI];
777  subProc[newBFaceI] = sourceProc[oldBFaceI];
778  subNewProc[newBFaceI] = sourceNewProc[oldBFaceI];
779  }
780  }
781 }
782 
783 
784 // Find cells on mesh whose faceID/procID match the neighbour cell/proc of
785 // domainMesh. Store the matching face.
786 void Foam::fvMeshDistribute::findCouples
787 (
788  const primitiveMesh& mesh,
789  const labelList& sourceFace,
790  const labelList& sourceProc,
791 
792  const label domain,
793  const primitiveMesh& domainMesh,
794  const labelList& domainFace,
795  const labelList& domainProc,
796 
797  labelList& masterCoupledFaces,
798  labelList& slaveCoupledFaces
799 )
800 {
801  // Store domain neighbour as map so we can easily look for pair
802  // with same face+proc.
803  HashTable<label, labelPair, labelPair::Hash<> > map(domainFace.size());
804 
805  forAll(domainFace, bFaceI)
806  {
807  map.insert(labelPair(domainFace[bFaceI], domainProc[bFaceI]), bFaceI);
808  }
809 
810 
811  // Try to match mesh data.
812 
813  masterCoupledFaces.setSize(domainFace.size());
814  slaveCoupledFaces.setSize(domainFace.size());
815  label coupledI = 0;
816 
817  forAll(sourceFace, bFaceI)
818  {
819  if (sourceProc[bFaceI] != -1)
820  {
821  labelPair myData(sourceFace[bFaceI], sourceProc[bFaceI]);
822 
823  HashTable<label, labelPair, labelPair::Hash<> >::const_iterator
824  iter = map.find(myData);
825 
826  if (iter != map.end())
827  {
828  label nbrBFaceI = iter();
829 
830  masterCoupledFaces[coupledI] = mesh.nInternalFaces() + bFaceI;
831  slaveCoupledFaces[coupledI] =
832  domainMesh.nInternalFaces()
833  + nbrBFaceI;
834 
835  coupledI++;
836  }
837  }
838  }
839 
840  masterCoupledFaces.setSize(coupledI);
841  slaveCoupledFaces.setSize(coupledI);
842 
843  if (debug)
844  {
845  Pout<< "findCouples : found " << coupledI
846  << " faces that will be stitched" << nl << endl;
847  }
848 }
849 
850 
851 // Map data on boundary faces to new mesh (resulting from adding two meshes)
852 Foam::labelList Foam::fvMeshDistribute::mapBoundaryData
853 (
854  const primitiveMesh& mesh, // mesh after adding
855  const mapAddedPolyMesh& map,
856  const labelList& boundaryData0, // mesh before adding
857  const label nInternalFaces1,
858  const labelList& boundaryData1 // added mesh
859 )
860 {
861  labelList newBoundaryData(mesh.nFaces() - mesh.nInternalFaces());
862 
863  forAll(boundaryData0, oldBFaceI)
864  {
865  label newFaceI = map.oldFaceMap()[oldBFaceI + map.nOldInternalFaces()];
866 
867  // Face still exists (is necessary?) and still boundary face
868  if (newFaceI >= 0 && newFaceI >= mesh.nInternalFaces())
869  {
870  newBoundaryData[newFaceI - mesh.nInternalFaces()] =
871  boundaryData0[oldBFaceI];
872  }
873  }
874 
875  forAll(boundaryData1, addedBFaceI)
876  {
877  label newFaceI = map.addedFaceMap()[addedBFaceI + nInternalFaces1];
878 
879  if (newFaceI >= 0 && newFaceI >= mesh.nInternalFaces())
880  {
881  newBoundaryData[newFaceI - mesh.nInternalFaces()] =
882  boundaryData1[addedBFaceI];
883  }
884  }
885 
886  return newBoundaryData;
887 }
888 
889 
890 // Remove cells. Add all exposed faces to patch oldInternalPatchI
891 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::doRemoveCells
892 (
893  const labelList& cellsToRemove,
894  const label oldInternalPatchI
895 )
896 {
897  // Mesh change engine
898  polyTopoChange meshMod(mesh_);
899 
900  // Cell removal topo engine. Do NOT synchronize parallel since
901  // we are doing a local cell removal.
902  removeCells cellRemover(mesh_, false);
903 
904  // Get all exposed faces
905  labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
906 
907  // Insert the topo changes
908  cellRemover.setRefinement
909  (
910  cellsToRemove,
911  exposedFaces,
912  labelList(exposedFaces.size(), oldInternalPatchI), // patch for exposed
913  // faces.
914  meshMod
915  );
916 
917  // Change the mesh. No inflation. Note: no parallel comms allowed.
918  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, false);
919 
920  // Update fields
921  mesh_.updateMesh(map);
922 
923  // Move mesh (since morphing does not do this)
924  if (map().hasMotionPoints())
925  {
926  mesh_.movePoints(map().preMotionPoints());
927  }
928 
929  return map;
930 }
931 
932 
933 // Delete and add processor patches. Changes mesh and returns per neighbour proc
934 // the processor patchID.
935 void Foam::fvMeshDistribute::addProcPatches
936 (
937  const labelList& neighbourNewProc, // processor that neighbour is on
938  labelList& procPatchID
939 )
940 {
941  // Now use the neighbourFace/Proc to repatch the mesh. These two lists
942  // contain for all current boundary faces the global patchID (for non-proc
943  // patch) or the processor.
944 
945  labelList procPatchSizes(Pstream::nProcs(), 0);
946 
947  forAll(neighbourNewProc, bFaceI)
948  {
949  if (neighbourNewProc[bFaceI] != -1)
950  {
951  procPatchSizes[neighbourNewProc[bFaceI]]++;
952  }
953  }
954 
955  // Per neighbour processor the label of the processor patch
956  procPatchID.setSize(Pstream::nProcs());
957 
958  forAll(procPatchSizes, procI)
959  {
960  if (procPatchSizes[procI] > 0)
961  {
962  const word patchName =
963  "procBoundary"
965  + "to"
966  + name(procI);
967 
968 
969  procPatchID[procI] = addProcPatch(patchName, procI);
970  addPatchFields<volScalarField>
971  (
972  processorFvPatchField<scalar>::typeName
973  );
974  addPatchFields<volVectorField>
975  (
976  processorFvPatchField<vector>::typeName
977  );
978  addPatchFields<volSphericalTensorField>
979  (
980  processorFvPatchField<sphericalTensor>::typeName
981  );
982  addPatchFields<volSymmTensorField>
983  (
984  processorFvPatchField<symmTensor>::typeName
985  );
986  addPatchFields<volTensorField>
987  (
988  processorFvPatchField<tensor>::typeName
989  );
990 
991  addPatchFields<surfaceScalarField>
992  (
993  processorFvPatchField<scalar>::typeName
994  );
995  addPatchFields<surfaceVectorField>
996  (
997  processorFvPatchField<vector>::typeName
998  );
999  addPatchFields<surfaceSphericalTensorField>
1000  (
1001  processorFvPatchField<sphericalTensor>::typeName
1002  );
1003  addPatchFields<surfaceSymmTensorField>
1004  (
1005  processorFvPatchField<symmTensor>::typeName
1006  );
1007  addPatchFields<surfaceTensorField>
1008  (
1009  processorFvPatchField<tensor>::typeName
1010  );
1011  }
1012  else
1013  {
1014  procPatchID[procI] = -1;
1015  }
1016  }
1017 }
1018 
1019 
1020 // Get boundary faces to be repatched. Is -1 or new patchID
1021 Foam::labelList Foam::fvMeshDistribute::getProcBoundaryPatch
1022 (
1023  const labelList& neighbourNewProc, // new processor per boundary face
1024  const labelList& procPatchID // patchID
1025 )
1026 {
1027  labelList patchIDs(neighbourNewProc);
1028 
1029  forAll(neighbourNewProc, bFaceI)
1030  {
1031  if (neighbourNewProc[bFaceI] != -1)
1032  {
1033  label nbrProc = neighbourNewProc[bFaceI];
1034 
1035  patchIDs[bFaceI] = procPatchID[nbrProc];
1036  }
1037  else
1038  {
1039  patchIDs[bFaceI] = -1;
1040  }
1041  }
1042  return patchIDs;
1043 }
1044 
1045 
1046 // Send mesh and coupling data.
1047 void Foam::fvMeshDistribute::sendMesh
1048 (
1049  const label domain,
1050  const fvMesh& mesh,
1051 
1052  const wordList& pointZoneNames,
1053  const wordList& faceZoneNames,
1054  const wordList& cellZoneNames,
1055 
1056  const labelList& sourceFace,
1057  const labelList& sourceProc,
1058  const labelList& sourceNewProc,
1059  OSstream& toDomain
1060 )
1061 {
1062  if (debug)
1063  {
1064  Pout<< "Sending to domain " << domain << nl
1065  << " nPoints:" << mesh.nPoints() << nl
1066  << " nFaces:" << mesh.nFaces() << nl
1067  << " nCells:" << mesh.nCells() << nl
1068  << " nPatches:" << mesh.boundaryMesh().size() << nl
1069  << endl;
1070  }
1071 
1072  // Assume sparse, possibly overlapping point zones. Get contents
1073  // in merged-zone indices.
1074  CompactListList_dev<label> zonePoints;
1075  {
1076  const pointZoneMesh& pointZones = mesh.pointZones();
1077 
1078  labelList rowSizes(pointZoneNames.size(), 0);
1079 
1080  forAll(pointZoneNames, nameI)
1081  {
1082  label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
1083 
1084  if (myZoneID != -1)
1085  {
1086  rowSizes[nameI] = pointZones[myZoneID].size();
1087  }
1088  }
1089  zonePoints.setSize(rowSizes);
1090 
1091  forAll(pointZoneNames, nameI)
1092  {
1093  label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
1094 
1095  if (myZoneID != -1)
1096  {
1097  zonePoints[nameI].assign(pointZones[myZoneID]);
1098  }
1099  }
1100  }
1101 
1102  // Assume sparse, possibly overlapping face zones
1103  CompactListList_dev<label> zoneFaces;
1104  CompactListList_dev<bool> zoneFaceFlip;
1105  {
1106  const faceZoneMesh& faceZones = mesh.faceZones();
1107 
1108  labelList rowSizes(faceZoneNames.size(), 0);
1109 
1110  forAll(faceZoneNames, nameI)
1111  {
1112  label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
1113 
1114  if (myZoneID != -1)
1115  {
1116  rowSizes[nameI] = faceZones[myZoneID].size();
1117  }
1118  }
1119 
1120  zoneFaces.setSize(rowSizes);
1121  zoneFaceFlip.setSize(rowSizes);
1122 
1123  forAll(faceZoneNames, nameI)
1124  {
1125  label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
1126 
1127  if (myZoneID != -1)
1128  {
1129  zoneFaces[nameI].assign(faceZones[myZoneID]);
1130  zoneFaceFlip[nameI].assign(faceZones[myZoneID].flipMap());
1131  }
1132  }
1133  }
1134 
1135  // Assume sparse, possibly overlapping cell zones
1136  CompactListList_dev<label> zoneCells;
1137  {
1138  const cellZoneMesh& cellZones = mesh.cellZones();
1139 
1140  labelList rowSizes(cellZoneNames.size(), 0);
1141 
1142  forAll(cellZoneNames, nameI)
1143  {
1144  label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
1145 
1146  if (myZoneID != -1)
1147  {
1148  rowSizes[nameI] = cellZones[myZoneID].size();
1149  }
1150  }
1151 
1152  zoneCells.setSize(rowSizes);
1153 
1154  forAll(cellZoneNames, nameI)
1155  {
1156  label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
1157 
1158  if (myZoneID != -1)
1159  {
1160  zoneCells[nameI].assign(cellZones[myZoneID]);
1161  }
1162  }
1163  }
1165  //labelList cellZoneID;
1166  //if (hasCellZones)
1167  //{
1168  // cellZoneID.setSize(mesh.nCells());
1169  // cellZoneID = -1;
1170  //
1171  // const cellZoneMesh& cellZones = mesh.cellZones();
1172  //
1173  // forAll(cellZones, zoneI)
1174  // {
1175  // UIndirectList<label>(cellZoneID, cellZones[zoneI]) = zoneI;
1176  // }
1177  //}
1178 
1179  // Send
1180  toDomain
1181  << mesh.points()
1182  << CompactListList_dev<label, face>(mesh.faces())
1183  << mesh.faceOwner()
1184  << mesh.faceNeighbour()
1185  << mesh.boundaryMesh()
1186 
1187  << zonePoints
1188  << zoneFaces
1189  << zoneFaceFlip
1190  << zoneCells
1191 
1192  << sourceFace
1193  << sourceProc
1194  << sourceNewProc;
1195 
1196 
1197  if (debug)
1198  {
1199  Pout<< "Started sending mesh to domain " << domain
1200  << endl;
1201  }
1202 }
1203 
1204 
1205 // Receive mesh. Opposite of sendMesh
1206 Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
1207 (
1208  const label domain,
1209  const wordList& pointZoneNames,
1210  const wordList& faceZoneNames,
1211  const wordList& cellZoneNames,
1212  const Time& runTime,
1213  labelList& domainSourceFace,
1214  labelList& domainSourceProc,
1215  labelList& domainSourceNewProc,
1216  ISstream& fromNbr
1217 )
1218 {
1219  pointField domainPoints(fromNbr);
1220  faceList domainFaces = CompactListList_dev<label, face>(fromNbr)();
1221  labelList domainAllOwner(fromNbr);
1222  labelList domainAllNeighbour(fromNbr);
1223  PtrList<entry> patchEntries(fromNbr);
1224 
1225  CompactListList_dev<label> zonePoints(fromNbr);
1226  CompactListList_dev<label> zoneFaces(fromNbr);
1227  CompactListList_dev<bool> zoneFaceFlip(fromNbr);
1228  CompactListList_dev<label> zoneCells(fromNbr);
1229 
1230  fromNbr
1231  >> domainSourceFace
1232  >> domainSourceProc
1233  >> domainSourceNewProc;
1234 
1235  // Construct fvMesh
1236  autoPtr<fvMesh> domainMeshPtr
1237  (
1238  new fvMesh
1239  (
1240  IOobject
1241  (
1243  runTime.timeName(),
1244  runTime,
1246  ),
1247  xferMove(domainPoints),
1248  xferMove(domainFaces),
1249  xferMove(domainAllOwner),
1250  xferMove(domainAllNeighbour),
1251  false // no parallel comms
1252  )
1253  );
1254  fvMesh& domainMesh = domainMeshPtr();
1255 
1256  List<polyPatch*> patches(patchEntries.size());
1257 
1258  forAll(patchEntries, patchI)
1259  {
1260  patches[patchI] = polyPatch::New
1261  (
1262  patchEntries[patchI].keyword(),
1263  patchEntries[patchI].dict(),
1264  patchI,
1265  domainMesh.boundaryMesh()
1266  ).ptr();
1267  }
1268  // Add patches; no parallel comms
1269  domainMesh.addFvPatches(patches, false);
1270 
1271  // Construct zones
1272  List<pointZone*> pZonePtrs(pointZoneNames.size());
1273  forAll(pZonePtrs, i)
1274  {
1275  pZonePtrs[i] = new pointZone
1276  (
1277  pointZoneNames[i],
1278  zonePoints[i],
1279  i,
1280  domainMesh.pointZones()
1281  );
1282  }
1283 
1284  List<faceZone*> fZonePtrs(faceZoneNames.size());
1285  forAll(fZonePtrs, i)
1286  {
1287  fZonePtrs[i] = new faceZone
1288  (
1289  faceZoneNames[i],
1290  zoneFaces[i],
1291  zoneFaceFlip[i],
1292  i,
1293  domainMesh.faceZones()
1294  );
1295  }
1296 
1297  List<cellZone*> cZonePtrs(cellZoneNames.size());
1298  forAll(cZonePtrs, i)
1299  {
1300  cZonePtrs[i] = new cellZone
1301  (
1302  cellZoneNames[i],
1303  zoneCells[i],
1304  i,
1305  domainMesh.cellZones()
1306  );
1307  }
1308  domainMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
1309 
1310  return domainMeshPtr;
1311 }
1312 
1313 
1314 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1315 
1316 // Construct from components
1317 Foam::fvMeshDistribute::fvMeshDistribute(fvMesh& mesh, const scalar mergeTol)
1318 :
1319  mesh_(mesh),
1320  mergeTol_(mergeTol)
1321 {}
1322 
1323 
1324 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1325 
1328  const labelList& distribution
1329 )
1330 {
1331  labelList nCells(Pstream::nProcs(), 0);
1332  forAll(distribution, cellI)
1333  {
1334  label newProc = distribution[cellI];
1335 
1336  if (newProc < 0 || newProc >= Pstream::nProcs())
1337  {
1338  FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
1339  << "Distribution should be in range 0.." << Pstream::nProcs()-1
1340  << endl
1341  << "At index " << cellI << " distribution:" << newProc
1342  << abort(FatalError);
1343  }
1344  nCells[newProc]++;
1345  }
1346  return nCells;
1347 }
1348 
1349 
1352  const labelList& distribution
1353 )
1354 {
1355  // Some checks on distribution
1356  if (distribution.size() != mesh_.nCells())
1357  {
1358  FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
1359  << "Size of distribution:"
1360  << distribution.size() << " mesh nCells:" << mesh_.nCells()
1361  << abort(FatalError);
1362  }
1363 
1364 
1365  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1366 
1367  // Check all processors have same non-proc patches in same order.
1368  if (patches.checkParallelSync(true))
1369  {
1370  FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
1371  << "This application requires all non-processor patches"
1372  << " to be present in the same order on all patches" << nl
1373  << "followed by the processor patches (which of course are unique)."
1374  << nl
1375  << "Local patches:" << mesh_.boundaryMesh().names()
1376  << abort(FatalError);
1377  }
1378 
1379  // Save some data for mapping later on
1380  const label nOldPoints(mesh_.nPoints());
1381  const label nOldFaces(mesh_.nFaces());
1382  const label nOldCells(mesh_.nCells());
1383  labelList oldPatchStarts(patches.size());
1384  labelList oldPatchNMeshPoints(patches.size());
1385  forAll(patches, patchI)
1386  {
1387  oldPatchStarts[patchI] = patches[patchI].start();
1388  oldPatchNMeshPoints[patchI] = patches[patchI].nPoints();
1389  }
1390 
1391 
1392 
1393  // Short circuit trivial case.
1394  if (!Pstream::parRun())
1395  {
1396  // Collect all maps and return
1398  (
1400  (
1401  mesh_,
1402 
1403  nOldPoints,
1404  nOldFaces,
1405  nOldCells,
1406  oldPatchStarts.xfer(),
1407  oldPatchNMeshPoints.xfer(),
1408 
1409  labelListList(1, identity(mesh_.nPoints())).xfer(),//subPointMap
1410  labelListList(1, identity(mesh_.nFaces())).xfer(), //subFaceMap
1411  labelListList(1, identity(mesh_.nCells())).xfer(), //subCellMap
1412  labelListList(1, identity(patches.size())).xfer(), //subPatchMap
1413 
1414  labelListList(1, identity(mesh_.nPoints())).xfer(),//pointMap
1415  labelListList(1, identity(mesh_.nFaces())).xfer(), //faceMap
1416  labelListList(1, identity(mesh_.nCells())).xfer(), //cellMap
1417  labelListList(1, identity(patches.size())).xfer() //patchMap
1418  )
1419  );
1420  }
1421 
1422 
1423  // Collect any zone names
1424  const wordList pointZoneNames(mergeWordList(mesh_.pointZones().names()));
1425  const wordList faceZoneNames(mergeWordList(mesh_.faceZones().names()));
1426  const wordList cellZoneNames(mergeWordList(mesh_.cellZones().names()));
1427 
1428 
1429 
1430  // Local environment of all boundary faces
1431  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1432 
1433  // A face is uniquely defined by
1434  // - proc
1435  // - local face no
1436  //
1437  // To glue the parts of meshes which can get sent from anywhere we
1438  // need to know on boundary faces what the above tuple on both sides is.
1439  // So we need to maintain:
1440  // - original face
1441  // - original processor id (= trivial)
1442  // For coupled boundaries (where the faces are 'duplicate') we take the
1443  // lowest numbered processor as the data to store.
1444  //
1445  // Additionally to create the procboundaries we need to know where the owner
1446  // cell on the other side now is: newNeighbourProc.
1447  //
1448 
1449  // physical boundary:
1450  // sourceProc = -1
1451  // sourceNewProc = -1
1452  // sourceFace = patchID
1453  // coupled boundary:
1454  // sourceProc = proc
1455  // sourceNewProc = distribution[cell on proc]
1456  // sourceFace = face
1457  labelList sourceFace;
1458  labelList sourceProc;
1459  labelList sourceNewProc;
1460  getNeighbourData(distribution, sourceFace, sourceProc, sourceNewProc);
1461 
1462 
1463  // Remove meshPhi. Since this would otherwise dissappear anyway
1464  // during topo changes and we have to guarantee that all the fields
1465  // can be sent.
1466  mesh_.clearOut();
1467  mesh_.resetMotion();
1468 
1469  // Get data to send. Make sure is synchronised
1470  const wordList volScalars(mesh_.names(volScalarField::typeName));
1471  checkEqualWordList("volScalarFields", volScalars);
1472  const wordList volVectors(mesh_.names(volVectorField::typeName));
1473  checkEqualWordList("volVectorFields", volVectors);
1474  const wordList volSphereTensors
1475  (
1477  );
1478  checkEqualWordList("volSphericalTensorFields", volSphereTensors);
1479  const wordList volSymmTensors(mesh_.names(volSymmTensorField::typeName));
1480  checkEqualWordList("volSymmTensorFields", volSymmTensors);
1481  const wordList volTensors(mesh_.names(volTensorField::typeName));
1482  checkEqualWordList("volTensorField", volTensors);
1483 
1484  const wordList surfScalars(mesh_.names(surfaceScalarField::typeName));
1485  checkEqualWordList("surfaceScalarFields", surfScalars);
1486  const wordList surfVectors(mesh_.names(surfaceVectorField::typeName));
1487  checkEqualWordList("surfaceVectorFields", surfVectors);
1488  const wordList surfSphereTensors
1489  (
1491  );
1492  checkEqualWordList("surfaceSphericalTensorFields", surfSphereTensors);
1493  const wordList surfSymmTensors
1494  (
1496  );
1497  checkEqualWordList("surfaceSymmTensorFields", surfSymmTensors);
1498  const wordList surfTensors(mesh_.names(surfaceTensorField::typeName));
1499  checkEqualWordList("surfaceTensorFields", surfTensors);
1500 
1501 
1502 
1503 
1504  // Find patch to temporarily put exposed and processor faces into.
1505  label oldInternalPatchI = findNonEmptyPatch();
1506 
1507 
1508 
1509  // Delete processor patches, starting from the back. Move all faces into
1510  // oldInternalPatchI.
1511  labelList repatchFaceMap;
1512  {
1513  autoPtr<mapPolyMesh> repatchMap = deleteProcPatches(oldInternalPatchI);
1514 
1515  // Store face map (only face ordering that changed)
1516  repatchFaceMap = repatchMap().faceMap();
1517 
1518  // Reorder all boundary face data (sourceProc, sourceFace etc.)
1519  labelList bFaceMap
1520  (
1522  (
1523  repatchMap().reverseFaceMap(),
1524  mesh_.nFaces() - mesh_.nInternalFaces(),
1525  mesh_.nInternalFaces()
1526  )
1527  - mesh_.nInternalFaces()
1528  );
1529 
1530  inplaceReorder(bFaceMap, sourceFace);
1531  inplaceReorder(bFaceMap, sourceProc);
1532  inplaceReorder(bFaceMap, sourceNewProc);
1533  }
1534 
1535 
1536 
1537  // Print a bit.
1538  if (debug)
1539  {
1540  Pout<< nl << "MESH WITH PROC PATCHES DELETED:" << endl;
1541  printMeshInfo(mesh_);
1542  printFieldInfo<volScalarField>(mesh_);
1543  printFieldInfo<volVectorField>(mesh_);
1544  printFieldInfo<volSphericalTensorField>(mesh_);
1545  printFieldInfo<volSymmTensorField>(mesh_);
1546  printFieldInfo<volTensorField>(mesh_);
1547  printFieldInfo<surfaceScalarField>(mesh_);
1548  printFieldInfo<surfaceVectorField>(mesh_);
1549  printFieldInfo<surfaceSphericalTensorField>(mesh_);
1550  printFieldInfo<surfaceSymmTensorField>(mesh_);
1551  printFieldInfo<surfaceTensorField>(mesh_);
1552  Pout<< nl << endl;
1553  }
1554 
1555 
1556 
1557  // Maps from subsetted mesh (that is sent) back to original maps
1558  labelListList subCellMap(Pstream::nProcs());
1559  labelListList subFaceMap(Pstream::nProcs());
1560  labelListList subPointMap(Pstream::nProcs());
1561  labelListList subPatchMap(Pstream::nProcs());
1562  // Maps from subsetted mesh to reconstructed mesh
1563  labelListList constructCellMap(Pstream::nProcs());
1564  labelListList constructFaceMap(Pstream::nProcs());
1565  labelListList constructPointMap(Pstream::nProcs());
1566  labelListList constructPatchMap(Pstream::nProcs());
1567 
1568 
1569 
1570 
1571  // Find out schedule
1572  // ~~~~~~~~~~~~~~~~~
1573 
1574  labelListList nSendCells(Pstream::nProcs());
1575  nSendCells[Pstream::myProcNo()] = countCells(distribution);
1576  Pstream::gatherList(nSendCells);
1577  Pstream::scatterList(nSendCells);
1578 
1579 
1580  // Allocate buffers
1582  forAll(sendStr, i)
1583  {
1584  sendStr.set(i, new OStringStream(IOstream::BINARY));
1585  }
1586  //PstreamBuffers pBufs(Pstream::nonBlocking);
1587 
1588 
1589  // What to send to neighbouring domains
1590  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1591 
1592  forAll(nSendCells[Pstream::myProcNo()], recvProc)
1593  {
1594  if
1595  (
1596  recvProc != Pstream::myProcNo()
1597  && nSendCells[Pstream::myProcNo()][recvProc] > 0
1598  )
1599  {
1600  // Send to recvProc
1601 
1602  if (debug)
1603  {
1604  Pout<< nl
1605  << "SUBSETTING FOR DOMAIN " << recvProc
1606  << " cells to send:"
1607  << nSendCells[Pstream::myProcNo()][recvProc]
1608  << nl << endl;
1609  }
1610 
1611  // Pstream for sending mesh and fields
1612  //OPstream str(Pstream::blocking, recvProc);
1613  //UOPstream str(recvProc, pBufs);
1614 
1615  // Mesh subsetting engine
1616  fvMeshSubset subsetter(mesh_);
1617 
1618  // Subset the cells of the current domain.
1619  subsetter.setLargeCellSubset
1620  (
1621  distribution,
1622  recvProc,
1623  oldInternalPatchI, // oldInternalFaces patch
1624  false // no parallel sync
1625  );
1626 
1627  subCellMap[recvProc] = subsetter.cellMap();
1628  subFaceMap[recvProc] = renumber
1629  (
1630  repatchFaceMap,
1631  subsetter.faceMap()
1632  );
1633  subPointMap[recvProc] = subsetter.pointMap();
1634  subPatchMap[recvProc] = subsetter.patchMap();
1635 
1636 
1637  // Subset the boundary fields (owner/neighbour/processor)
1638  labelList procSourceFace;
1639  labelList procSourceProc;
1640  labelList procSourceNewProc;
1641 
1642  subsetBoundaryData
1643  (
1644  subsetter.subMesh(),
1645  subsetter.faceMap(), // from subMesh to mesh
1646  subsetter.cellMap(), // ,, ,,
1647 
1648  distribution, // old mesh distribution
1649  mesh_.faceOwner(), // old owner
1650  mesh_.faceNeighbour(),
1651  mesh_.nInternalFaces(),
1652 
1653  sourceFace,
1654  sourceProc,
1655  sourceNewProc,
1656 
1657  procSourceFace,
1658  procSourceProc,
1659  procSourceNewProc
1660  );
1661 
1662 
1663 
1664  // Send to neighbour
1665  sendMesh
1666  (
1667  recvProc,
1668  subsetter.subMesh(),
1669 
1670  pointZoneNames,
1671  faceZoneNames,
1672  cellZoneNames,
1673 
1674  procSourceFace,
1675  procSourceProc,
1676  procSourceNewProc,
1677  sendStr[recvProc]
1678  );
1679  sendFields<volScalarField>
1680  (
1681  recvProc,
1682  volScalars,
1683  subsetter,
1684  sendStr[recvProc]
1685  );
1686  sendFields<volVectorField>
1687  (
1688  recvProc,
1689  volVectors,
1690  subsetter,
1691  sendStr[recvProc]
1692  );
1693  sendFields<volSphericalTensorField>
1694  (
1695  recvProc,
1696  volSphereTensors,
1697  subsetter,
1698  sendStr[recvProc]
1699  );
1700  sendFields<volSymmTensorField>
1701  (
1702  recvProc,
1703  volSymmTensors,
1704  subsetter,
1705  sendStr[recvProc]
1706  );
1707  sendFields<volTensorField>
1708  (
1709  recvProc,
1710  volTensors,
1711  subsetter,
1712  sendStr[recvProc]
1713  );
1714 
1715  sendFields<surfaceScalarField>
1716  (
1717  recvProc,
1718  surfScalars,
1719  subsetter,
1720  sendStr[recvProc]
1721  );
1722  sendFields<surfaceVectorField>
1723  (
1724  recvProc,
1725  surfVectors,
1726  subsetter,
1727  sendStr[recvProc]
1728  );
1729  sendFields<surfaceSphericalTensorField>
1730  (
1731  recvProc,
1732  surfSphereTensors,
1733  subsetter,
1734  sendStr[recvProc]
1735  );
1736  sendFields<surfaceSymmTensorField>
1737  (
1738  recvProc,
1739  surfSymmTensors,
1740  subsetter,
1741  sendStr[recvProc]
1742  );
1743  sendFields<surfaceTensorField>
1744  (
1745  recvProc,
1746  surfTensors,
1747  subsetter,
1748  sendStr[recvProc]
1749  );
1750  }
1751  }
1752 
1753 
1754  // Start sending&receiving from buffers
1755  //pBufs.finishedSends();
1756 
1757  // get the data.
1759  {
1760  List<List<char> > sendBufs(sendStr.size());
1761  forAll(sendStr, procI)
1762  {
1763  string contents = sendStr[procI].str();
1764  const char* ptr = contents.data();
1765 
1766  sendBufs[procI].setSize(contents.size());
1767  forAll(sendBufs[procI], i)
1768  {
1769  sendBufs[procI][i] = *ptr++;
1770  }
1771  // Clear OStringStream
1772  sendStr.set(procI, NULL);
1773  }
1774 
1775  // Transfer sendBufs into recvBufs
1776  List<List<char> > recvBufs(Pstream::nProcs());
1777  labelListList sizes(Pstream::nProcs());
1778  exchange<List<char>, char>(sendBufs, recvBufs, sizes);
1779 
1780  forAll(recvStr, procI)
1781  {
1782  string contents(recvBufs[procI].begin(), recvBufs[procI].size());
1783  recvStr.set
1784  (
1785  procI,
1786  new IStringStream(contents, IOstream::BINARY)
1787  );
1788  }
1789  }
1790 
1791 
1792  // Subset the part that stays
1793  // ~~~~~~~~~~~~~~~~~~~~~~~~~~
1794 
1795  {
1796  // Save old mesh maps before changing mesh
1797  const labelList oldFaceOwner(mesh_.faceOwner());
1798  const labelList oldFaceNeighbour(mesh_.faceNeighbour());
1799  const label oldInternalFaces = mesh_.nInternalFaces();
1800 
1801  // Remove cells.
1802  autoPtr<mapPolyMesh> subMap
1803  (
1804  doRemoveCells
1805  (
1806  select(false, distribution, Pstream::myProcNo()),
1807  oldInternalPatchI
1808  )
1809  );
1810 
1811  // Addressing from subsetted mesh
1812  subCellMap[Pstream::myProcNo()] = subMap().cellMap();
1813  subFaceMap[Pstream::myProcNo()] = renumber
1814  (
1815  repatchFaceMap,
1816  subMap().faceMap()
1817  );
1818  subPointMap[Pstream::myProcNo()] = subMap().pointMap();
1819  subPatchMap[Pstream::myProcNo()] = identity(patches.size());
1820 
1821  // Initialize all addressing into current mesh
1822  constructCellMap[Pstream::myProcNo()] = identity(mesh_.nCells());
1823  constructFaceMap[Pstream::myProcNo()] = identity(mesh_.nFaces());
1824  constructPointMap[Pstream::myProcNo()] = identity(mesh_.nPoints());
1825  constructPatchMap[Pstream::myProcNo()] = identity(patches.size());
1826 
1827  // Subset the mesh data: neighbourCell/neighbourProc
1828  // fields
1829  labelList domainSourceFace;
1830  labelList domainSourceProc;
1831  labelList domainSourceNewProc;
1832 
1833  subsetBoundaryData
1834  (
1835  mesh_, // new mesh
1836  subMap().faceMap(), // from new to original mesh
1837  subMap().cellMap(),
1838 
1839  distribution, // distribution before subsetting
1840  oldFaceOwner, // owner before subsetting
1841  oldFaceNeighbour, // neighbour ,,
1842  oldInternalFaces, // nInternalFaces ,,
1843 
1844  sourceFace,
1845  sourceProc,
1846  sourceNewProc,
1847 
1848  domainSourceFace,
1849  domainSourceProc,
1850  domainSourceNewProc
1851  );
1852 
1853  sourceFace.transfer(domainSourceFace);
1854  sourceProc.transfer(domainSourceProc);
1855  sourceNewProc.transfer(domainSourceNewProc);
1856  }
1857 
1858 
1859  // Print a bit.
1860  if (debug)
1861  {
1862  Pout<< nl << "STARTING MESH:" << endl;
1863  printMeshInfo(mesh_);
1864  printFieldInfo<volScalarField>(mesh_);
1865  printFieldInfo<volVectorField>(mesh_);
1866  printFieldInfo<volSphericalTensorField>(mesh_);
1867  printFieldInfo<volSymmTensorField>(mesh_);
1868  printFieldInfo<volTensorField>(mesh_);
1869  printFieldInfo<surfaceScalarField>(mesh_);
1870  printFieldInfo<surfaceVectorField>(mesh_);
1871  printFieldInfo<surfaceSphericalTensorField>(mesh_);
1872  printFieldInfo<surfaceSymmTensorField>(mesh_);
1873  printFieldInfo<surfaceTensorField>(mesh_);
1874  Pout<< nl << endl;
1875  }
1876 
1877 
1878 
1879  // Receive and add what was sent
1880  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1881 
1882  forAll(nSendCells, sendProc)
1883  {
1884  // Did processor sendProc send anything to me?
1885  if
1886  (
1887  sendProc != Pstream::myProcNo()
1888  && nSendCells[sendProc][Pstream::myProcNo()] > 0
1889  )
1890  {
1891  if (debug)
1892  {
1893  Pout<< nl
1894  << "RECEIVING FROM DOMAIN " << sendProc
1895  << " cells to receive:"
1896  << nSendCells[sendProc][Pstream::myProcNo()]
1897  << nl << endl;
1898  }
1899 
1900 
1901  // Pstream for receiving mesh and fields
1902  //UIPstream str(sendProc, pBufs);
1903 
1904 
1905  // Receive from sendProc
1906  labelList domainSourceFace;
1907  labelList domainSourceProc;
1908  labelList domainSourceNewProc;
1909 
1910  autoPtr<fvMesh> domainMeshPtr;
1921 
1922  // Opposite of sendMesh
1923  {
1924  domainMeshPtr = receiveMesh
1925  (
1926  sendProc,
1927  pointZoneNames,
1928  faceZoneNames,
1929  cellZoneNames,
1930 
1931  const_cast<Time&>(mesh_.time()),
1932  domainSourceFace,
1933  domainSourceProc,
1934  domainSourceNewProc,
1935  recvStr[sendProc]
1936  );
1937  fvMesh& domainMesh = domainMeshPtr();
1938 
1939  // Receive fields. Read as single dictionary because
1940  // of problems reading consecutive fields from single stream.
1941  dictionary fieldDicts(recvStr[sendProc]);
1942 
1943  receiveFields<volScalarField>
1944  (
1945  sendProc,
1946  volScalars,
1947  domainMesh,
1948  vsf,
1949  fieldDicts.subDict(volScalarField::typeName)
1950  );
1951  receiveFields<volVectorField>
1952  (
1953  sendProc,
1954  volVectors,
1955  domainMesh,
1956  vvf,
1957  fieldDicts.subDict(volVectorField::typeName)
1958  );
1959  receiveFields<volSphericalTensorField>
1960  (
1961  sendProc,
1962  volSphereTensors,
1963  domainMesh,
1964  vsptf,
1966  );
1967  receiveFields<volSymmTensorField>
1968  (
1969  sendProc,
1970  volSymmTensors,
1971  domainMesh,
1972  vsytf,
1974  );
1975  receiveFields<volTensorField>
1976  (
1977  sendProc,
1978  volTensors,
1979  domainMesh,
1980  vtf,
1981  fieldDicts.subDict(volTensorField::typeName)
1982  );
1983 
1984  receiveFields<surfaceScalarField>
1985  (
1986  sendProc,
1987  surfScalars,
1988  domainMesh,
1989  ssf,
1991  );
1992  receiveFields<surfaceVectorField>
1993  (
1994  sendProc,
1995  surfVectors,
1996  domainMesh,
1997  svf,
1999  );
2000  receiveFields<surfaceSphericalTensorField>
2001  (
2002  sendProc,
2003  surfSphereTensors,
2004  domainMesh,
2005  ssptf,
2007  );
2008  receiveFields<surfaceSymmTensorField>
2009  (
2010  sendProc,
2011  surfSymmTensors,
2012  domainMesh,
2013  ssytf,
2015  );
2016  receiveFields<surfaceTensorField>
2017  (
2018  sendProc,
2019  surfTensors,
2020  domainMesh,
2021  stf,
2023  );
2024  }
2025  const fvMesh& domainMesh = domainMeshPtr();
2026 
2027 
2028  constructCellMap[sendProc] = identity(domainMesh.nCells());
2029  constructFaceMap[sendProc] = identity(domainMesh.nFaces());
2030  constructPointMap[sendProc] = identity(domainMesh.nPoints());
2031  constructPatchMap[sendProc] =
2032  identity(domainMesh.boundaryMesh().size());
2033 
2034 
2035  // Print a bit.
2036  if (debug)
2037  {
2038  Pout<< nl << "RECEIVED MESH FROM:" << sendProc << endl;
2039  printMeshInfo(domainMesh);
2040  printFieldInfo<volScalarField>(domainMesh);
2041  printFieldInfo<volVectorField>(domainMesh);
2042  printFieldInfo<volSphericalTensorField>(domainMesh);
2043  printFieldInfo<volSymmTensorField>(domainMesh);
2044  printFieldInfo<volTensorField>(domainMesh);
2045  printFieldInfo<surfaceScalarField>(domainMesh);
2046  printFieldInfo<surfaceVectorField>(domainMesh);
2047  printFieldInfo<surfaceSphericalTensorField>(domainMesh);
2048  printFieldInfo<surfaceSymmTensorField>(domainMesh);
2049  printFieldInfo<surfaceTensorField>(domainMesh);
2050  }
2051 
2052 
2053  // Now this mesh we received (from sendProc) needs to be merged
2054  // with the current mesh. On the current mesh we have for all
2055  // boundaryfaces the original face and processor. See if we can
2056  // match these up to the received domainSourceFace and
2057  // domainSourceProc.
2058  labelList masterCoupledFaces;
2059  labelList slaveCoupledFaces;
2060  findCouples
2061  (
2062  mesh_,
2063 
2064  sourceFace,
2065  sourceProc,
2066 
2067  sendProc,
2068  domainMesh,
2069  domainSourceFace,
2070  domainSourceProc,
2071 
2072  masterCoupledFaces,
2073  slaveCoupledFaces
2074  );
2075 
2076  // Generate additional coupling info (points, edges) from
2077  // faces-that-match
2078  faceCoupleInfo couples
2079  (
2080  mesh_,
2081  masterCoupledFaces,
2082  domainMesh,
2083  slaveCoupledFaces,
2084  mergeTol_, // merge tolerance
2085  true, // faces align
2086  true, // couples are ordered already
2087  false
2088  );
2089 
2090 
2091  // Add domainMesh to mesh
2092  // ~~~~~~~~~~~~~~~~~~~~~~
2093 
2095  (
2096  mesh_,
2097  domainMesh,
2098  couples,
2099  false // no parallel comms
2100  );
2101 
2102  // Update mesh data: sourceFace,sourceProc for added
2103  // mesh.
2104 
2105  sourceFace =
2106  mapBoundaryData
2107  (
2108  mesh_,
2109  map(),
2110  sourceFace,
2111  domainMesh.nInternalFaces(),
2112  domainSourceFace
2113  );
2114  sourceProc =
2115  mapBoundaryData
2116  (
2117  mesh_,
2118  map(),
2119  sourceProc,
2120  domainMesh.nInternalFaces(),
2121  domainSourceProc
2122  );
2123  sourceNewProc =
2124  mapBoundaryData
2125  (
2126  mesh_,
2127  map(),
2128  sourceNewProc,
2129  domainMesh.nInternalFaces(),
2130  domainSourceNewProc
2131  );
2132 
2133  // Update all addressing so xxProcAddressing points to correct item
2134  // in masterMesh.
2135  const labelList& oldCellMap = map().oldCellMap();
2136  const labelList& oldFaceMap = map().oldFaceMap();
2137  const labelList& oldPointMap = map().oldPointMap();
2138  const labelList& oldPatchMap = map().oldPatchMap();
2139 
2140  forAll(constructPatchMap, procI)
2141  {
2142  if (procI != sendProc && constructPatchMap[procI].size())
2143  {
2144  // Processor already in mesh (either myProcNo or received)
2145  inplaceRenumber(oldCellMap, constructCellMap[procI]);
2146  inplaceRenumber(oldFaceMap, constructFaceMap[procI]);
2147  inplaceRenumber(oldPointMap, constructPointMap[procI]);
2148  inplaceRenumber(oldPatchMap, constructPatchMap[procI]);
2149  }
2150  }
2151 
2152  // Added processor
2153  inplaceRenumber(map().addedCellMap(), constructCellMap[sendProc]);
2154  inplaceRenumber(map().addedFaceMap(), constructFaceMap[sendProc]);
2155  inplaceRenumber(map().addedPointMap(), constructPointMap[sendProc]);
2156  inplaceRenumber(map().addedPatchMap(), constructPatchMap[sendProc]);
2157 
2158  if (debug)
2159  {
2160  Pout<< nl << "MERGED MESH FROM:" << sendProc << endl;
2161  printMeshInfo(mesh_);
2162  printFieldInfo<volScalarField>(mesh_);
2163  printFieldInfo<volVectorField>(mesh_);
2164  printFieldInfo<volSphericalTensorField>(mesh_);
2165  printFieldInfo<volSymmTensorField>(mesh_);
2166  printFieldInfo<volTensorField>(mesh_);
2167  printFieldInfo<surfaceScalarField>(mesh_);
2168  printFieldInfo<surfaceVectorField>(mesh_);
2169  printFieldInfo<surfaceSphericalTensorField>(mesh_);
2170  printFieldInfo<surfaceSymmTensorField>(mesh_);
2171  printFieldInfo<surfaceTensorField>(mesh_);
2172  Pout<< nl << endl;
2173  }
2174  }
2175  }
2176 
2177 
2178  // Print a bit.
2179  if (debug)
2180  {
2181  Pout<< nl << "REDISTRIBUTED MESH:" << endl;
2182  printMeshInfo(mesh_);
2183  printFieldInfo<volScalarField>(mesh_);
2184  printFieldInfo<volVectorField>(mesh_);
2185  printFieldInfo<volSphericalTensorField>(mesh_);
2186  printFieldInfo<volSymmTensorField>(mesh_);
2187  printFieldInfo<volTensorField>(mesh_);
2188  printFieldInfo<surfaceScalarField>(mesh_);
2189  printFieldInfo<surfaceVectorField>(mesh_);
2190  printFieldInfo<surfaceSphericalTensorField>(mesh_);
2191  printFieldInfo<surfaceSymmTensorField>(mesh_);
2192  printFieldInfo<surfaceTensorField>(mesh_);
2193  Pout<< nl << endl;
2194  }
2195 
2196 
2197 
2198  // Add processorPatches
2199  // ~~~~~~~~~~~~~~~~~~~~
2200 
2201  // Per neighbour processor the patchID to it (or -1).
2202  labelList procPatchID;
2203 
2204  // Add processor patches.
2205  addProcPatches(sourceNewProc, procPatchID);
2206 
2207  // Put faces into correct patch. Note that we now have proper
2208  // processorPolyPatches again so repatching will take care of coupled face
2209  // ordering.
2210 
2211  // Get boundary faces to be repatched. Is -1 or new patchID
2212  labelList newPatchID
2213  (
2214  getProcBoundaryPatch
2215  (
2216  sourceNewProc,
2217  procPatchID
2218  )
2219  );
2220 
2221  // Change patches. Since this might change ordering of coupled faces
2222  // we also need to adapt our constructMaps.
2223  repatch(newPatchID, constructFaceMap);
2224 
2225  // See if any geometrically shared points need to be merged. Note: does
2226  // parallel comms.
2227  mergeSharedPoints(constructPointMap);
2228 
2229  // Bit of hack: processorFvPatchField does not get reset since created
2230  // from nothing so explicitly reset.
2231  initPatchFields<volScalarField, processorFvPatchField<scalar> >
2232  (
2234  );
2235  initPatchFields<volVectorField, processorFvPatchField<vector> >
2236  (
2238  );
2239  initPatchFields
2240  <
2243  >
2244  (
2246  );
2247  initPatchFields<volSymmTensorField, processorFvPatchField<symmTensor> >
2248  (
2250  );
2251  initPatchFields<volTensorField, processorFvPatchField<tensor> >
2252  (
2254  );
2255  initPatchFields<surfaceScalarField, processorFvsPatchField<scalar> >
2256  (
2258  );
2259  initPatchFields<surfaceVectorField, processorFvsPatchField<vector> >
2260  (
2262  );
2263  initPatchFields
2264  <
2267  >
2268  (
2270  );
2271  initPatchFields
2272  <
2275  >
2276  (
2278  );
2279  initPatchFields<surfaceTensorField, processorFvsPatchField<tensor> >
2280  (
2282  );
2283 
2284 
2285  mesh_.setInstance(mesh_.time().timeName());
2286 
2287 
2288  // Print a bit
2289  if (debug)
2290  {
2291  Pout<< nl << "FINAL MESH:" << endl;
2292  printMeshInfo(mesh_);
2293  printFieldInfo<volScalarField>(mesh_);
2294  printFieldInfo<volVectorField>(mesh_);
2295  printFieldInfo<volSphericalTensorField>(mesh_);
2296  printFieldInfo<volSymmTensorField>(mesh_);
2297  printFieldInfo<volTensorField>(mesh_);
2298  printFieldInfo<surfaceScalarField>(mesh_);
2299  printFieldInfo<surfaceVectorField>(mesh_);
2300  printFieldInfo<surfaceSphericalTensorField>(mesh_);
2301  printFieldInfo<surfaceSymmTensorField>(mesh_);
2302  printFieldInfo<surfaceTensorField>(mesh_);
2303  Pout<< nl << endl;
2304  }
2305 
2306  // Collect all maps and return
2308  (
2310  (
2311  mesh_,
2312 
2313  nOldPoints,
2314  nOldFaces,
2315  nOldCells,
2316  oldPatchStarts,
2317  oldPatchNMeshPoints,
2318 
2319  subPointMap,
2320  subFaceMap,
2321  subCellMap,
2322  subPatchMap,
2323 
2324  constructPointMap,
2325  constructFaceMap,
2326  constructCellMap,
2327  constructPatchMap,
2328  true // reuse storage
2329  )
2330  );
2331 }
2332 
2333 
2334 // ************************ vim: set sw=4 sts=4 et: ************************ //