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