FreeFOAM The Cross-Platform CFD Toolkit
fvMeshSubset.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 Description
25  Post-processing mesh subset tool. Given the original mesh and the
26  list of selected cells, it creates the mesh consisting only of the
27  desired cells, with the mapping list for points, faces, and cells.
28 
29  MJ 23/03/05 on coupled faces change the patch of the face to the
30  oldInternalFaces patch.
31 
32 \*---------------------------------------------------------------------------*/
33 
34 #include "fvMeshSubset.H"
35 #include <OpenFOAM/boolList.H>
36 #include <OpenFOAM/Pstream.H>
40 
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 
43 namespace Foam
44 {
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 bool Foam::fvMeshSubset::checkCellSubset() const
49 {
50  if (fvMeshSubsetPtr_.empty())
51  {
52  FatalErrorIn("bool fvMeshSubset::checkCellSubset() const")
53  << "Mesh subset not set. Please set the cell map using "
54  << "void setCellSubset(const labelHashSet& cellsToSubset)" << endl
55  << "before attempting to access subset data"
56  << abort(FatalError);
57 
58  return false;
59  }
60  else
61  {
62  return true;
63  }
64 }
65 
66 
67 void Foam::fvMeshSubset::markPoints
68 (
69  const labelList& curPoints,
70  Map<label>& pointMap
71 )
72 {
73  forAll (curPoints, pointI)
74  {
75  // Note: insert will only insert if not yet there.
76  pointMap.insert(curPoints[pointI], 0);
77  }
78 }
79 
80 
81 void Foam::fvMeshSubset::markPoints
82 (
83  const labelList& curPoints,
84  labelList& pointMap
85 )
86 {
87  forAll (curPoints, pointI)
88  {
89  pointMap[curPoints[pointI]] = 0;
90  }
91 }
92 
93 
94 // Synchronize nCellsUsingFace on both sides of coupled patches. Marks
95 // faces that become 'uncoupled' with 3.
96 void Foam::fvMeshSubset::doCoupledPatches
97 (
98  const bool syncPar,
99  labelList& nCellsUsingFace
100 ) const
101 {
102  const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
103 
104  label nUncoupled = 0;
105 
106  if (syncPar && Pstream::parRun())
107  {
108  // Send face usage across processor patches
109  forAll (oldPatches, oldPatchI)
110  {
111  const polyPatch& pp = oldPatches[oldPatchI];
112 
113  if (isA<processorPolyPatch>(pp))
114  {
115  const processorPolyPatch& procPatch =
116  refCast<const processorPolyPatch>(pp);
117 
118  OPstream toNeighbour
119  (
121  procPatch.neighbProcNo()
122  );
123 
124  toNeighbour
125  << SubList<label>(nCellsUsingFace, pp.size(), pp.start());
126  }
127  }
128 
129 
130  // Receive face usage count and check for faces that become uncoupled.
131  forAll (oldPatches, oldPatchI)
132  {
133  const polyPatch& pp = oldPatches[oldPatchI];
134 
135  if (isA<processorPolyPatch>(pp))
136  {
137  const processorPolyPatch& procPatch =
138  refCast<const processorPolyPatch>(pp);
139 
140  IPstream fromNeighbour
141  (
143  procPatch.neighbProcNo()
144  );
145 
146  labelList nbrCellsUsingFace(fromNeighbour);
147 
148  // Combine with this side.
149 
150  forAll (pp, i)
151  {
152  if
153  (
154  nCellsUsingFace[pp.start()+i] == 1
155  && nbrCellsUsingFace[i] == 0
156  )
157  {
158  // Face's neighbour is no longer there. Mark face off
159  // as coupled
160  nCellsUsingFace[pp.start()+i] = 3;
161  nUncoupled++;
162  }
163  }
164  }
165  }
166  }
167 
168  // Do same for cyclics.
169  forAll (oldPatches, oldPatchI)
170  {
171  const polyPatch& pp = oldPatches[oldPatchI];
172 
173  if (isA<cyclicPolyPatch>(pp))
174  {
175  const cyclicPolyPatch& cycPatch =
176  refCast<const cyclicPolyPatch>(pp);
177 
178  forAll (cycPatch, i)
179  {
180  label thisFaceI = cycPatch.start() + i;
181  label otherFaceI = cycPatch.transformGlobalFace(thisFaceI);
182 
183  if
184  (
185  nCellsUsingFace[thisFaceI] == 1
186  && nCellsUsingFace[otherFaceI] == 0
187  )
188  {
189  nCellsUsingFace[thisFaceI] = 3;
190  nUncoupled++;
191  }
192  }
193  }
194  }
195 
196  if (syncPar)
197  {
198  reduce(nUncoupled, sumOp<label>());
199  }
200 
201  if (nUncoupled > 0)
202  {
203  Info<< "Uncoupled " << nUncoupled << " faces on coupled patches. "
204  << "(processorPolyPatch, cyclicPolyPatch)" << nl
205  << "You might need to run couplePatches to restore the patch face"
206  << " ordering." << nl << endl;
207  }
208 }
209 
210 
212 (
213  const label nElems,
214  const labelList& selectedElements,
215  const labelList& subsetMap
216 )
217 {
218  // Mark selected elements.
219  boolList selected(nElems, false);
220  forAll(selectedElements, i)
221  {
222  selected[selectedElements[i]] = true;
223  }
224 
225  // Count subset of selected elements
226  label n = 0;
227  forAll(subsetMap, i)
228  {
229  if (selected[subsetMap[i]])
230  {
231  n++;
232  }
233  }
234 
235  // Collect selected elements
236  labelList subsettedElements(n);
237  n = 0;
238 
239  forAll(subsetMap, i)
240  {
241  if (selected[subsetMap[i]])
242  {
243  subsettedElements[n++] = i;
244  }
245  }
246 
247  return subsettedElements;
248 }
249 
250 
251 void Foam::fvMeshSubset::subsetZones()
252 {
253  // Keep all zones, even if zero size.
254 
255  const pointZoneMesh& pointZones = baseMesh().pointZones();
256 
257  // PointZones
258  List<pointZone*> pZonePtrs(pointZones.size());
259 
260  forAll(pointZones, i)
261  {
262  const pointZone& pz = pointZones[i];
263 
264  pZonePtrs[i] = new pointZone
265  (
266  pz.name(),
267  subset(baseMesh().nPoints(), pz, pointMap()),
268  i,
269  fvMeshSubsetPtr_().pointZones()
270  );
271  }
272 
273 
274  // FaceZones
275 
276  const faceZoneMesh& faceZones = baseMesh().faceZones();
277 
278 
279  // Do we need to remove zones where the side we're interested in
280  // no longer exists? Guess not.
281  List<faceZone*> fZonePtrs(faceZones.size());
282 
283  forAll(faceZones, i)
284  {
285  const faceZone& fz = faceZones[i];
286 
287  // Expand faceZone to full mesh
288  // +1 : part of faceZone, flipped
289  // -1 : ,, , unflipped
290  // 0 : not part of faceZone
291  labelList zone(baseMesh().nFaces(), 0);
292  forAll(fz, j)
293  {
294  if (fz.flipMap()[j])
295  {
296  zone[fz[j]] = 1;
297  }
298  else
299  {
300  zone[fz[j]] = -1;
301  }
302  }
303 
304  // Select faces
305  label nSub = 0;
306  forAll(faceMap(), j)
307  {
308  if (zone[faceMap()[j]] != 0)
309  {
310  nSub++;
311  }
312  }
313  labelList subAddressing(nSub);
314  boolList subFlipStatus(nSub);
315  nSub = 0;
316  forAll(faceMap(), subFaceI)
317  {
318  label meshFaceI = faceMap()[subFaceI];
319  if (zone[meshFaceI] != 0)
320  {
321  subAddressing[nSub] = subFaceI;
322  label subOwner = subMesh().faceOwner()[subFaceI];
323  label baseOwner = baseMesh().faceOwner()[meshFaceI];
324  // If subowner is the same cell as the base keep the flip status
325  bool sameOwner = (cellMap()[subOwner] == baseOwner);
326  bool flip = (zone[meshFaceI] == 1);
327  subFlipStatus[nSub] = (sameOwner == flip);
328 
329  nSub++;
330  }
331  }
332 
333  fZonePtrs[i] = new faceZone
334  (
335  fz.name(),
336  subAddressing,
337  subFlipStatus,
338  i,
339  fvMeshSubsetPtr_().faceZones()
340  );
341  }
342 
343 
344  const cellZoneMesh& cellZones = baseMesh().cellZones();
345 
346  List<cellZone*> cZonePtrs(cellZones.size());
347 
348  forAll(cellZones, i)
349  {
350  const cellZone& cz = cellZones[i];
351 
352  cZonePtrs[i] = new cellZone
353  (
354  cz.name(),
355  subset(baseMesh().nCells(), cz, cellMap()),
356  i,
357  fvMeshSubsetPtr_().cellZones()
358  );
359  }
360 
361 
362  // Add the zones
363  fvMeshSubsetPtr_().addZones(pZonePtrs, fZonePtrs, cZonePtrs);
364 }
365 
366 
367 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
368 
369 // Construct from components
370 Foam::fvMeshSubset::fvMeshSubset(const fvMesh& baseMesh)
371 :
372  baseMesh_(baseMesh),
373  fvMeshSubsetPtr_(NULL),
374  pointMap_(0),
375  faceMap_(0),
376  cellMap_(0),
377  patchMap_(0)
378 {}
379 
380 
381 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
382 
384 (
385  const labelHashSet& globalCellMap,
386  const label patchID
387 )
388 {
389  // Initial check on patches before doing anything time consuming.
390  const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
391  const cellList& oldCells = baseMesh().cells();
392  const faceList& oldFaces = baseMesh().faces();
393  const pointField& oldPoints = baseMesh().points();
394  const labelList& oldOwner = baseMesh().faceOwner();
395  const labelList& oldNeighbour = baseMesh().faceNeighbour();
396 
397  label wantedPatchID = patchID;
398 
399  if (wantedPatchID == -1)
400  {
401  // No explicit patch specified. Put in oldInternalFaces patch.
402  // Check if patch with this name already exists.
403  wantedPatchID = oldPatches.findPatchID("oldInternalFaces");
404  }
405  else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.size())
406  {
408  (
409  "fvMeshSubset::setCellSubset(const labelHashSet&"
410  ", const label patchID)"
411  ) << "Non-existing patch index " << wantedPatchID << endl
412  << "Should be between 0 and " << oldPatches.size()-1
413  << abort(FatalError);
414  }
415 
416 
417  cellMap_ = globalCellMap.toc();
418 
419  // Sort the cell map in the ascending order
420  sort(cellMap_);
421 
422  // Approximate sizing parameters for face and point lists
423  const label avgNFacesPerCell = 6;
424  const label avgNPointsPerFace = 4;
425 
426 
427  label nCellsInSet = cellMap_.size();
428 
429  // Mark all used faces
430 
431  Map<label> facesToSubset(avgNFacesPerCell*nCellsInSet);
432 
433  forAll (cellMap_, cellI)
434  {
435  // Mark all faces from the cell
436  const labelList& curFaces = oldCells[cellMap_[cellI]];
437 
438  forAll (curFaces, faceI)
439  {
440  if (!facesToSubset.found(curFaces[faceI]))
441  {
442  facesToSubset.insert(curFaces[faceI], 1);
443  }
444  else
445  {
446  facesToSubset[curFaces[faceI]]++;
447  }
448  }
449  }
450 
451  // Mark all used points and make a global-to-local face map
452  Map<label> globalFaceMap(facesToSubset.size());
453 
454  // Make a global-to-local point map
455  Map<label> globalPointMap(avgNPointsPerFace*facesToSubset.size());
456 
457  // This is done in two goes, so that the boundary faces are last
458  // in the list. Because of this, I need to create the face map
459  // along the way rather than just grab the table of contents.
460  labelList facesToc = facesToSubset.toc();
461  sort(facesToc);
462  faceMap_.setSize(facesToc.size());
463 
464  // 1. Get all faces that will be internal to the submesh.
465  forAll (facesToc, faceI)
466  {
467  if (facesToSubset[facesToc[faceI]] == 2)
468  {
469  // Mark face and increment number of points in set
470  faceMap_[globalFaceMap.size()] = facesToc[faceI];
471  globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
472 
473  // Mark all points from the face
474  markPoints(oldFaces[facesToc[faceI]], globalPointMap);
475  }
476  }
477 
478  // These are all the internal faces in the mesh.
479  label nInternalFaces = globalFaceMap.size();
480 
481 
482  // Where to insert old internal faces.
483  label oldPatchStart = labelMax;
484  if (wantedPatchID != -1)
485  {
486  oldPatchStart = oldPatches[wantedPatchID].start();
487  }
488 
489 
490  label faceI = 0;
491 
492  // 2. Boundary faces up to where we want to insert old internal faces
493  for (; faceI< facesToc.size(); faceI++)
494  {
495  if (facesToc[faceI] >= oldPatchStart)
496  {
497  break;
498  }
499  if
500  (
501  !baseMesh().isInternalFace(facesToc[faceI])
502  && facesToSubset[facesToc[faceI]] == 1
503  )
504  {
505  // Mark face and increment number of points in set
506  faceMap_[globalFaceMap.size()] = facesToc[faceI];
507  globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
508 
509  // Mark all points from the face
510  markPoints(oldFaces[facesToc[faceI]], globalPointMap);
511  }
512  }
513 
514  // 3. old internal faces
515  forAll(facesToc, intFaceI)
516  {
517  if
518  (
519  baseMesh().isInternalFace(facesToc[intFaceI])
520  && facesToSubset[facesToc[intFaceI]] == 1
521  )
522  {
523  // Mark face and increment number of points in set
524  faceMap_[globalFaceMap.size()] = facesToc[intFaceI];
525  globalFaceMap.insert(facesToc[intFaceI], globalFaceMap.size());
526 
527  // Mark all points from the face
528  markPoints(oldFaces[facesToc[intFaceI]], globalPointMap);
529  }
530  }
531 
532  // 4. Remaining boundary faces
533  for (; faceI< facesToc.size(); faceI++)
534  {
535  if
536  (
537  !baseMesh().isInternalFace(facesToc[faceI])
538  && facesToSubset[facesToc[faceI]] == 1
539  )
540  {
541  // Mark face and increment number of points in set
542  faceMap_[globalFaceMap.size()] = facesToc[faceI];
543  globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
544 
545  // Mark all points from the face
546  markPoints(oldFaces[facesToc[faceI]], globalPointMap);
547  }
548  }
549 
550 
551 
552  // Grab the points map
553  pointMap_ = globalPointMap.toc();
554  sort(pointMap_);
555 
556  forAll (pointMap_, pointI)
557  {
558  globalPointMap[pointMap_[pointI]] = pointI;
559  }
560 
561  Pout << "Number of cells in new mesh: " << nCellsInSet << endl;
562  Pout << "Number of faces in new mesh: " << globalFaceMap.size() << endl;
563  Pout << "Number of points in new mesh: " << globalPointMap.size() << endl;
564 
565  // Make a new mesh
566  pointField newPoints(globalPointMap.size());
567 
568  label nNewPoints = 0;
569 
570  forAll (pointMap_, pointI)
571  {
572  newPoints[nNewPoints] = oldPoints[pointMap_[pointI]];
573  nNewPoints++;
574  }
575 
576  faceList newFaces(globalFaceMap.size());
577 
578  label nNewFaces = 0;
579 
580  // Make internal faces
581  for (label faceI = 0; faceI < nInternalFaces; faceI++)
582  {
583  const face& oldF = oldFaces[faceMap_[faceI]];
584 
585  face newF(oldF.size());
586 
587  forAll (newF, i)
588  {
589  newF[i] = globalPointMap[oldF[i]];
590  }
591 
592  newFaces[nNewFaces] = newF;
593  nNewFaces++;
594  }
595 
596  // Make boundary faces
597 
598  label nbSize = oldPatches.size();
599  label oldInternalPatchID = -1;
600 
601  if (wantedPatchID == -1)
602  {
603  // Create 'oldInternalFaces' patch at the end
604  // and put all exposed internal faces in there.
605  oldInternalPatchID = nbSize;
606  nbSize++;
607 
608  }
609  else
610  {
611  oldInternalPatchID = wantedPatchID;
612  }
613 
614 
615  // Grad size and start of each patch on the fly. Because of the
616  // structure of the underlying mesh, the patches will appear in the
617  // ascending order
618  labelList boundaryPatchSizes(nbSize, 0);
619 
620  // Assign boundary faces. Visited in order of faceMap_.
621  for (label faceI = nInternalFaces; faceI < faceMap_.size(); faceI++)
622  {
623  label oldFaceI = faceMap_[faceI];
624 
625  face oldF = oldFaces[oldFaceI];
626 
627  // Turn the faces as necessary to point outwards
628  if (baseMesh().isInternalFace(oldFaceI))
629  {
630  // Internal face. Possibly turned the wrong way round
631  if
632  (
633  !globalCellMap.found(oldOwner[oldFaceI])
634  && globalCellMap.found(oldNeighbour[oldFaceI])
635  )
636  {
637  oldF = oldFaces[oldFaceI].reverseFace();
638  }
639 
640  // Update count for patch
641  boundaryPatchSizes[oldInternalPatchID]++;
642  }
643  else
644  {
645  // Boundary face. Increment the appropriate patch
646  label patchOfFace = oldPatches.whichPatch(oldFaceI);
647 
648  // Update count for patch
649  boundaryPatchSizes[patchOfFace]++;
650  }
651 
652  face newF(oldF.size());
653 
654  forAll (newF, i)
655  {
656  newF[i] = globalPointMap[oldF[i]];
657  }
658 
659  newFaces[nNewFaces] = newF;
660  nNewFaces++;
661  }
662 
663 
664 
665  // Create cells
666  cellList newCells(nCellsInSet);
667 
668  label nNewCells = 0;
669 
670  forAll (cellMap_, cellI)
671  {
672  const labelList& oldC = oldCells[cellMap_[cellI]];
673 
674  labelList newC(oldC.size());
675 
676  forAll (newC, i)
677  {
678  newC[i] = globalFaceMap[oldC[i]];
679  }
680 
681  newCells[nNewCells] = cell(newC);
682  nNewCells++;
683  }
684 
685 
686  // Delete any old one
687  fvMeshSubsetPtr_.clear();
688  // Make a new mesh
689  fvMeshSubsetPtr_.reset
690  (
691  new fvMesh
692  (
693  IOobject
694  (
695  baseMesh().name() + "SubSet",
696  baseMesh().time().timeName(),
697  baseMesh().time(),
700  ),
701  xferMove(newPoints),
702  xferMove(newFaces),
703  xferMove(newCells)
704  )
705  );
706 
707 
708  // Add old patches
709  List<polyPatch*> newBoundary(nbSize);
710  patchMap_.setSize(nbSize);
711  label nNewPatches = 0;
712  label patchStart = nInternalFaces;
713 
714  forAll(oldPatches, patchI)
715  {
716  if (boundaryPatchSizes[patchI] > 0)
717  {
718  // Patch still exists. Add it
719  newBoundary[nNewPatches] = oldPatches[patchI].clone
720  (
721  fvMeshSubsetPtr_().boundaryMesh(),
722  nNewPatches,
723  boundaryPatchSizes[patchI],
724  patchStart
725  ).ptr();
726 
727  patchStart += boundaryPatchSizes[patchI];
728  patchMap_[nNewPatches] = patchI;
729  nNewPatches++;
730  }
731  }
732 
733  if (wantedPatchID == -1)
734  {
735  // Newly created patch so is at end. Check if any faces in it.
736  if (boundaryPatchSizes[oldInternalPatchID] > 0)
737  {
738  newBoundary[nNewPatches] = new emptyPolyPatch
739  (
740  "oldInternalFaces",
741  boundaryPatchSizes[oldInternalPatchID],
742  patchStart,
743  nNewPatches,
744  fvMeshSubsetPtr_().boundaryMesh()
745  );
746 
747  // The index for the first patch is -1 as it originates from
748  // the internal faces
749  patchMap_[nNewPatches] = -1;
750  nNewPatches++;
751  }
752  }
753 
754  // Reset the patch lists
755  newBoundary.setSize(nNewPatches);
756  patchMap_.setSize(nNewPatches);
757 
758  // Add the fvPatches
759  fvMeshSubsetPtr_().addFvPatches(newBoundary);
760 
761  // Subset and add any zones
762  subsetZones();
763 }
764 
765 
767 (
768  const labelList& region,
769  const label currentRegion,
770  const label patchID,
771  const bool syncPar
772 )
773 {
774  const cellList& oldCells = baseMesh().cells();
775  const faceList& oldFaces = baseMesh().faces();
776  const pointField& oldPoints = baseMesh().points();
777  const labelList& oldOwner = baseMesh().faceOwner();
778  const labelList& oldNeighbour = baseMesh().faceNeighbour();
779  const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
780  const label oldNInternalFaces = baseMesh().nInternalFaces();
781 
782  // Initial checks
783 
784  if (region.size() != oldCells.size())
785  {
787  (
788  "fvMeshSubset::setCellSubset(const labelList&"
789  ", const label, const label, const bool)"
790  ) << "Size of region " << region.size()
791  << " is not equal to number of cells in mesh " << oldCells.size()
792  << abort(FatalError);
793  }
794 
795 
796  label wantedPatchID = patchID;
797 
798  if (wantedPatchID == -1)
799  {
800  // No explicit patch specified. Put in oldInternalFaces patch.
801  // Check if patch with this name already exists.
802  wantedPatchID = oldPatches.findPatchID("oldInternalFaces");
803  }
804  else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.size())
805  {
807  (
808  "fvMeshSubset::setCellSubset(const labelList&"
809  ", const label, const label, const bool)"
810  ) << "Non-existing patch index " << wantedPatchID << endl
811  << "Should be between 0 and " << oldPatches.size()-1
812  << abort(FatalError);
813  }
814 
815 
816  // Get the cells for the current region.
817  cellMap_.setSize(oldCells.size());
818  label nCellsInSet = 0;
819 
820  forAll (region, oldCellI)
821  {
822  if (region[oldCellI] == currentRegion)
823  {
824  cellMap_[nCellsInSet++] = oldCellI;
825  }
826  }
827  cellMap_.setSize(nCellsInSet);
828 
829 
830  // Mark all used faces. Count number of cells using them
831  // 0: face not used anymore
832  // 1: face used by one cell, face becomes/stays boundary face
833  // 2: face still used and remains internal face
834  // 3: face coupled and used by one cell only (so should become normal,
835  // non-coupled patch face)
836  //
837  // Note that this is not really nessecary - but means we can size things
838  // correctly. Also makes handling coupled faces much easier.
839 
840  labelList nCellsUsingFace(oldFaces.size(), 0);
841 
842  label nFacesInSet = 0;
843  forAll (oldFaces, oldFaceI)
844  {
845  bool faceUsed = false;
846 
847  if (region[oldOwner[oldFaceI]] == currentRegion)
848  {
849  nCellsUsingFace[oldFaceI]++;
850  faceUsed = true;
851  }
852 
853  if
854  (
855  baseMesh().isInternalFace(oldFaceI)
856  && (region[oldNeighbour[oldFaceI]] == currentRegion)
857  )
858  {
859  nCellsUsingFace[oldFaceI]++;
860  faceUsed = true;
861  }
862 
863  if (faceUsed)
864  {
865  nFacesInSet++;
866  }
867  }
868  faceMap_.setSize(nFacesInSet);
869 
870  // Handle coupled faces. Modifies patch faces to be uncoupled to 3.
871  doCoupledPatches(syncPar, nCellsUsingFace);
872 
873 
874  // See which patch to use for exposed internal faces.
875  label oldInternalPatchID = 0;
876 
877  // Insert faces before which patch
878  label nextPatchID = oldPatches.size();
879 
880  // old to new patches
881  labelList globalPatchMap(oldPatches.size());
882 
883  // New patch size
884  label nbSize = oldPatches.size();
885 
886  if (wantedPatchID == -1)
887  {
888  // Create 'oldInternalFaces' patch at the end (or before
889  // processorPatches)
890  // and put all exposed internal faces in there.
891 
892  forAll(oldPatches, patchI)
893  {
894  if (isA<processorPolyPatch>(oldPatches[patchI]))
895  {
896  nextPatchID = patchI;
897  break;
898  }
899  oldInternalPatchID++;
900  }
901 
902  nbSize++;
903 
904  // adapt old to new patches for inserted patch
905  for (label oldPatchI = 0; oldPatchI < nextPatchID; oldPatchI++)
906  {
907  globalPatchMap[oldPatchI] = oldPatchI;
908  }
909  for
910  (
911  label oldPatchI = nextPatchID;
912  oldPatchI < oldPatches.size();
913  oldPatchI++
914  )
915  {
916  globalPatchMap[oldPatchI] = oldPatchI+1;
917  }
918  }
919  else
920  {
921  oldInternalPatchID = wantedPatchID;
922  nextPatchID = wantedPatchID+1;
923 
924  // old to new patches
925  globalPatchMap = identity(oldPatches.size());
926  }
927 
928  labelList boundaryPatchSizes(nbSize, 0);
929 
930 
931  // Make a global-to-local point map
932  labelList globalPointMap(oldPoints.size(), -1);
933 
934  labelList globalFaceMap(oldFaces.size(), -1);
935  label faceI = 0;
936 
937  // 1. Pick up all preserved internal faces.
938  for (label oldFaceI = 0; oldFaceI < oldNInternalFaces; oldFaceI++)
939  {
940  if (nCellsUsingFace[oldFaceI] == 2)
941  {
942  globalFaceMap[oldFaceI] = faceI;
943  faceMap_[faceI++] = oldFaceI;
944 
945  // Mark all points from the face
946  markPoints(oldFaces[oldFaceI], globalPointMap);
947  }
948  }
949 
950  // These are all the internal faces in the mesh.
951  label nInternalFaces = faceI;
952 
953  // 2. Boundary faces up to where we want to insert old internal faces
954  for
955  (
956  label oldPatchI = 0;
957  oldPatchI < oldPatches.size()
958  && oldPatchI < nextPatchID;
959  oldPatchI++
960  )
961  {
962  const polyPatch& oldPatch = oldPatches[oldPatchI];
963 
964  label oldFaceI = oldPatch.start();
965 
966  forAll (oldPatch, i)
967  {
968  if (nCellsUsingFace[oldFaceI] == 1)
969  {
970  // Boundary face is kept.
971 
972  // Mark face and increment number of points in set
973  globalFaceMap[oldFaceI] = faceI;
974  faceMap_[faceI++] = oldFaceI;
975 
976  // Mark all points from the face
977  markPoints(oldFaces[oldFaceI], globalPointMap);
978 
979  // Increment number of patch faces
980  boundaryPatchSizes[globalPatchMap[oldPatchI]]++;
981  }
982  oldFaceI++;
983  }
984  }
985 
986  // 3a. old internal faces that have become exposed.
987  for (label oldFaceI = 0; oldFaceI < oldNInternalFaces; oldFaceI++)
988  {
989  if (nCellsUsingFace[oldFaceI] == 1)
990  {
991  globalFaceMap[oldFaceI] = faceI;
992  faceMap_[faceI++] = oldFaceI;
993 
994  // Mark all points from the face
995  markPoints(oldFaces[oldFaceI], globalPointMap);
996 
997  // Increment number of patch faces
998  boundaryPatchSizes[oldInternalPatchID]++;
999  }
1000  }
1001 
1002  // 3b. coupled patch faces that have become uncoupled.
1003  for
1004  (
1005  label oldFaceI = oldNInternalFaces;
1006  oldFaceI < oldFaces.size();
1007  oldFaceI++
1008  )
1009  {
1010  if (nCellsUsingFace[oldFaceI] == 3)
1011  {
1012  globalFaceMap[oldFaceI] = faceI;
1013  faceMap_[faceI++] = oldFaceI;
1014 
1015  // Mark all points from the face
1016  markPoints(oldFaces[oldFaceI], globalPointMap);
1017 
1018  // Increment number of patch faces
1019  boundaryPatchSizes[oldInternalPatchID]++;
1020  }
1021  }
1022 
1023  // 4. Remaining boundary faces
1024  for
1025  (
1026  label oldPatchI = nextPatchID;
1027  oldPatchI < oldPatches.size();
1028  oldPatchI++
1029  )
1030  {
1031  const polyPatch& oldPatch = oldPatches[oldPatchI];
1032 
1033  label oldFaceI = oldPatch.start();
1034 
1035  forAll (oldPatch, i)
1036  {
1037  if (nCellsUsingFace[oldFaceI] == 1)
1038  {
1039  // Boundary face is kept.
1040 
1041  // Mark face and increment number of points in set
1042  globalFaceMap[oldFaceI] = faceI;
1043  faceMap_[faceI++] = oldFaceI;
1044 
1045  // Mark all points from the face
1046  markPoints(oldFaces[oldFaceI], globalPointMap);
1047 
1048  // Increment number of patch faces
1049  boundaryPatchSizes[globalPatchMap[oldPatchI]]++;
1050  }
1051  oldFaceI++;
1052  }
1053  }
1054 
1055  if (faceI != nFacesInSet)
1056  {
1057  FatalErrorIn
1058  (
1059  "fvMeshSubset::setCellSubset(const labelList&"
1060  ", const label, const label, const bool)"
1061  ) << "Problem" << abort(FatalError);
1062  }
1063 
1064 
1065  // Grab the points map
1066  label nPointsInSet = 0;
1067 
1068  forAll(globalPointMap, pointI)
1069  {
1070  if (globalPointMap[pointI] != -1)
1071  {
1072  nPointsInSet++;
1073  }
1074  }
1075  pointMap_.setSize(nPointsInSet);
1076 
1077  nPointsInSet = 0;
1078 
1079  forAll(globalPointMap, pointI)
1080  {
1081  if (globalPointMap[pointI] != -1)
1082  {
1083  pointMap_[nPointsInSet] = pointI;
1084  globalPointMap[pointI] = nPointsInSet;
1085  nPointsInSet++;
1086  }
1087  }
1088 
1089  //Pout<< "Number of cells in new mesh : " << cellMap_.size() << endl;
1090  //Pout<< "Number of faces in new mesh : " << faceMap_.size() << endl;
1091  //Pout<< "Number of points in new mesh: " << pointMap_.size() << endl;
1092 
1093  // Make a new mesh
1094  pointField newPoints(pointMap_.size());
1095 
1096  label nNewPoints = 0;
1097 
1098  forAll (pointMap_, pointI)
1099  {
1100  newPoints[nNewPoints] = oldPoints[pointMap_[pointI]];
1101  nNewPoints++;
1102  }
1103 
1104  faceList newFaces(faceMap_.size());
1105 
1106  label nNewFaces = 0;
1107 
1108  // Make internal faces
1109  for (label faceI = 0; faceI < nInternalFaces; faceI++)
1110  {
1111  const face& oldF = oldFaces[faceMap_[faceI]];
1112 
1113  face newF(oldF.size());
1114 
1115  forAll (newF, i)
1116  {
1117  newF[i] = globalPointMap[oldF[i]];
1118  }
1119 
1120  newFaces[nNewFaces] = newF;
1121  nNewFaces++;
1122  }
1123 
1124 
1125  // Make boundary faces. (different from internal since might need to be
1126  // flipped)
1127  for (label faceI = nInternalFaces; faceI < faceMap_.size(); faceI++)
1128  {
1129  label oldFaceI = faceMap_[faceI];
1130 
1131  face oldF = oldFaces[oldFaceI];
1132 
1133  // Turn the faces as necessary to point outwards
1134  if (baseMesh().isInternalFace(oldFaceI))
1135  {
1136  // Was internal face. Possibly turned the wrong way round
1137  if
1138  (
1139  region[oldOwner[oldFaceI]] != currentRegion
1140  && region[oldNeighbour[oldFaceI]] == currentRegion
1141  )
1142  {
1143  oldF = oldFaces[oldFaceI].reverseFace();
1144  }
1145  }
1146 
1147  // Relabel vertices of the (possibly turned) face.
1148  face newF(oldF.size());
1149 
1150  forAll (newF, i)
1151  {
1152  newF[i] = globalPointMap[oldF[i]];
1153  }
1154 
1155  newFaces[nNewFaces] = newF;
1156  nNewFaces++;
1157  }
1158 
1159 
1160 
1161  // Create cells
1162  cellList newCells(nCellsInSet);
1163 
1164  label nNewCells = 0;
1165 
1166  forAll (cellMap_, cellI)
1167  {
1168  const labelList& oldC = oldCells[cellMap_[cellI]];
1169 
1170  labelList newC(oldC.size());
1171 
1172  forAll (newC, i)
1173  {
1174  newC[i] = globalFaceMap[oldC[i]];
1175  }
1176 
1177  newCells[nNewCells] = cell(newC);
1178  nNewCells++;
1179  }
1180 
1181 
1182  // Make a new mesh
1183  // Note that mesh gets registered with same name as original mesh. This is
1184  // not proper but cannot be avoided since otherwise surfaceInterpolation
1185  // cannot find its fvSchemes (it will try to read e.g.
1186  // system/region0SubSet/fvSchemes)
1187  fvMeshSubsetPtr_.reset
1188  (
1189  new fvMesh
1190  (
1191  IOobject
1192  (
1193  baseMesh().name(),
1194  baseMesh().time().timeName(),
1195  baseMesh().time(),
1198  ),
1199  xferMove(newPoints),
1200  xferMove(newFaces),
1201  xferMove(newCells),
1202  syncPar // parallel synchronisation
1203  )
1204  );
1205 
1206  // Add old patches
1207  List<polyPatch*> newBoundary(nbSize);
1208  patchMap_.setSize(nbSize);
1209  label nNewPatches = 0;
1210  label patchStart = nInternalFaces;
1211 
1212 
1213  // For parallel: only remove patch if none of the processors has it.
1214  // This only gets done for patches before the one being inserted
1215  // (so patches < nextPatchID)
1216 
1217  // Get sum of patch sizes. Zero if patch can be deleted.
1218  labelList globalPatchSizes(boundaryPatchSizes);
1219  globalPatchSizes.setSize(nextPatchID);
1220 
1221  if (syncPar && Pstream::parRun())
1222  {
1223  // Get patch names (up to nextPatchID)
1225  patchNames[Pstream::myProcNo()] = oldPatches.names();
1226  patchNames[Pstream::myProcNo()].setSize(nextPatchID);
1227  Pstream::gatherList(patchNames);
1228  Pstream::scatterList(patchNames);
1229 
1230  // Get patch sizes (up to nextPatchID).
1231  // Note that up to nextPatchID the globalPatchMap is an identity so
1232  // no need to index through that.
1233  Pstream::listCombineGather(globalPatchSizes, plusEqOp<label>());
1234  Pstream::listCombineScatter(globalPatchSizes);
1235 
1236  // Now all processors have all the patchnames.
1237  // Decide: if all processors have the same patch names and size is zero
1238  // everywhere remove the patch.
1239  bool samePatches = true;
1240 
1241  for (label procI = 1; procI < patchNames.size(); procI++)
1242  {
1243  if (patchNames[procI] != patchNames[0])
1244  {
1245  samePatches = false;
1246  break;
1247  }
1248  }
1249 
1250  if (!samePatches)
1251  {
1252  // Patchnames not sync on all processors so disable removal of
1253  // zero sized patches.
1254  globalPatchSizes = labelMax;
1255  }
1256  }
1257 
1258 
1259  // Old patches
1260 
1261  for
1262  (
1263  label oldPatchI = 0;
1264  oldPatchI < oldPatches.size()
1265  && oldPatchI < nextPatchID;
1266  oldPatchI++
1267  )
1268  {
1269  label newSize = boundaryPatchSizes[globalPatchMap[oldPatchI]];
1270 
1271  // Clone (even if 0 size)
1272  newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
1273  (
1274  fvMeshSubsetPtr_().boundaryMesh(),
1275  nNewPatches,
1276  newSize,
1277  patchStart
1278  ).ptr();
1279 
1280  patchStart += newSize;
1281  patchMap_[nNewPatches] = oldPatchI; // compact patchMap
1282  nNewPatches++;
1283  }
1284 
1285  // Inserted patch
1286 
1287  if (wantedPatchID == -1)
1288  {
1289  label oldInternalSize = boundaryPatchSizes[oldInternalPatchID];
1290 
1291  if (syncPar)
1292  {
1293  reduce(oldInternalSize, sumOp<label>());
1294  }
1295 
1296  // Newly created patch so is at end. Check if any faces in it.
1297  if (oldInternalSize > 0)
1298  {
1299  newBoundary[nNewPatches] = new emptyPolyPatch
1300  (
1301  "oldInternalFaces",
1302  boundaryPatchSizes[oldInternalPatchID],
1303  patchStart,
1304  nNewPatches,
1305  fvMeshSubsetPtr_().boundaryMesh()
1306  );
1307 
1308  //Pout<< " oldInternalFaces : "
1309  // << boundaryPatchSizes[oldInternalPatchID] << endl;
1310 
1311  // The index for the first patch is -1 as it originates from
1312  // the internal faces
1313  patchStart += boundaryPatchSizes[oldInternalPatchID];
1314  patchMap_[nNewPatches] = -1;
1315  nNewPatches++;
1316  }
1317  }
1318 
1319  // Old patches
1320 
1321  for
1322  (
1323  label oldPatchI = nextPatchID;
1324  oldPatchI < oldPatches.size();
1325  oldPatchI++
1326  )
1327  {
1328  label newSize = boundaryPatchSizes[globalPatchMap[oldPatchI]];
1329 
1330  // Patch still exists. Add it
1331  newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
1332  (
1333  fvMeshSubsetPtr_().boundaryMesh(),
1334  nNewPatches,
1335  newSize,
1336  patchStart
1337  ).ptr();
1338 
1339  //Pout<< " " << oldPatches[oldPatchI].name() << " : "
1340  // << newSize << endl;
1341 
1342  patchStart += newSize;
1343  patchMap_[nNewPatches] = oldPatchI; // compact patchMap
1344  nNewPatches++;
1345  }
1346 
1347 
1348  // Reset the patch lists
1349  newBoundary.setSize(nNewPatches);
1350  patchMap_.setSize(nNewPatches);
1351 
1352 
1353  // Add the fvPatches
1354  fvMeshSubsetPtr_().addFvPatches(newBoundary, syncPar);
1355 
1356  // Subset and add any zones
1357  subsetZones();
1358 }
1359 
1360 
1363  const labelHashSet& globalCellMap,
1364  const label patchID,
1365  const bool syncPar
1366 )
1367 {
1368  labelList region(baseMesh().nCells(), 0);
1369 
1370  forAllConstIter (labelHashSet, globalCellMap, iter)
1371  {
1372  region[iter.key()] = 1;
1373  }
1374  setLargeCellSubset(region, 1, patchID, syncPar);
1375 }
1376 
1377 
1379 {
1380  checkCellSubset();
1381 
1382  return fvMeshSubsetPtr_();
1383 }
1384 
1385 
1387 {
1388  checkCellSubset();
1389 
1390  return fvMeshSubsetPtr_();
1391 }
1392 
1393 
1395 {
1396  checkCellSubset();
1397 
1398  return pointMap_;
1399 }
1400 
1401 
1403 {
1404  checkCellSubset();
1405 
1406  return faceMap_;
1407 }
1408 
1409 
1411 {
1412  checkCellSubset();
1413 
1414  return cellMap_;
1415 }
1416 
1417 
1419 {
1420  checkCellSubset();
1421 
1422  return patchMap_;
1423 }
1424 
1425 
1426 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1427 
1428 } // End namespace Foam
1429 
1430 // ************************ vim: set sw=4 sts=4 et: ************************ //