1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright held by original author
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
26 Post-processing mesh subset tool. Given the original mesh and the
27 list of selected cells, it creates the mesh consisting only of the
28 desired cells, with the mapping list for points, faces, and cells.
30 MJ 23/03/05 on coupled faces change the patch of the face to the
31 oldInternalFaces patch.
33 \*---------------------------------------------------------------------------*/
35 #include "fvMeshSubset.H"
38 #include "emptyPolyPatch.H"
39 #include "demandDrivenData.H"
40 #include "cyclicPolyPatch.H"
42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 bool Foam::fvMeshSubset::checkCellSubset() const
51 if (!fvMeshSubsetPtr_)
53 FatalErrorIn("bool fvMeshSubset::checkCellSubset() const")
54 << "Mesh subset not set. Please set the cell map using "
55 << "void setCellSubset(const labelHashSet& cellsToSubset)" << endl
56 << "before attempting to access subset data"
68 void Foam::fvMeshSubset::markPoints
70 const labelList& curPoints,
74 forAll (curPoints, pointI)
76 // Note: insert will only insert if not yet there.
77 pointMap.insert(curPoints[pointI], 0);
82 void Foam::fvMeshSubset::markPoints
84 const labelList& curPoints,
88 forAll (curPoints, pointI)
90 pointMap[curPoints[pointI]] = 0;
95 // Synchronize nCellsUsingFace on both sides of coupled patches. Marks
96 // faces that become 'uncoupled' with 3.
97 void Foam::fvMeshSubset::doCoupledPatches
100 labelList& nCellsUsingFace
103 const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
105 label nUncoupled = 0;
107 if (syncPar && Pstream::parRun())
109 // Send face usage across processor patches
110 forAll (oldPatches, oldPatchI)
112 const polyPatch& pp = oldPatches[oldPatchI];
114 if (isA<processorPolyPatch>(pp))
116 const processorPolyPatch& procPatch =
117 refCast<const processorPolyPatch>(pp);
122 procPatch.neighbProcNo()
126 << SubList<label>(nCellsUsingFace, pp.size(), pp.start());
131 // Receive face usage count and check for faces that become uncoupled.
132 forAll (oldPatches, oldPatchI)
134 const polyPatch& pp = oldPatches[oldPatchI];
136 if (isA<processorPolyPatch>(pp))
138 const processorPolyPatch& procPatch =
139 refCast<const processorPolyPatch>(pp);
141 IPstream fromNeighbour
144 procPatch.neighbProcNo()
147 labelList nbrCellsUsingFace(fromNeighbour);
149 // Combine with this side.
155 nCellsUsingFace[pp.start()+i] == 1
156 && nbrCellsUsingFace[i] == 0
159 // Face's neighbour is no longer there. Mark face off
161 nCellsUsingFace[pp.start()+i] = 3;
169 // Do same for cyclics.
170 forAll (oldPatches, oldPatchI)
172 const polyPatch& pp = oldPatches[oldPatchI];
174 if (isA<cyclicPolyPatch>(pp))
176 const cyclicPolyPatch& cycPatch =
177 refCast<const cyclicPolyPatch>(pp);
181 label thisFaceI = cycPatch.start() + i;
182 label otherFaceI = cycPatch.transformGlobalFace(thisFaceI);
186 nCellsUsingFace[thisFaceI] == 1
187 && nCellsUsingFace[otherFaceI] == 0
190 nCellsUsingFace[thisFaceI] = 3;
199 reduce(nUncoupled, sumOp<label>());
204 Info<< "Uncoupled " << nUncoupled << " faces on coupled patches. "
205 << "(processorPolyPatch, cyclicPolyPatch)" << nl
206 << "You might need to run couplePatches to restore the patch face"
207 << " ordering." << nl << endl;
212 labelList Foam::fvMeshSubset::subset
215 const labelList& selectedElements,
216 const labelList& subsetMap
219 // Mark selected elements.
220 boolList selected(nElems, false);
221 forAll (selectedElements, i)
223 selected[selectedElements[i]] = true;
226 // Count subset of selected elements
228 forAll (subsetMap, i)
230 if (selected[subsetMap[i]])
236 // Collect selected elements
237 labelList subsettedElements(n);
242 if (selected[subsetMap[i]])
244 subsettedElements[n++] = i;
248 return subsettedElements;
252 void Foam::fvMeshSubset::subsetZones()
254 // Keep all zones, even if zero size.
256 const pointZoneMesh& pointZones = baseMesh().pointZones();
259 List<pointZone*> pZonePtrs(pointZones.size());
261 forAll(pointZones, i)
263 const pointZone& pz = pointZones[i];
265 pZonePtrs[i] = new pointZone
268 subset(baseMesh().allPoints().size(), pz, pointMap()),
270 fvMeshSubsetPtr_->pointZones()
277 const faceZoneMesh& faceZones = baseMesh().faceZones();
280 // Do we need to remove zones where the side we're interested in
281 // no longer exists? Guess not.
282 List<faceZone*> fZonePtrs(faceZones.size());
286 const faceZone& fz = faceZones[i];
288 // Create list of mesh faces part of the new zone
289 labelList subAddressing
293 baseMesh().allFaces().size(),
299 // Flipmap for all mesh faces
300 boolList fullFlipStatus(baseMesh().allFaces().size(), false);
303 fullFlipStatus[fz[j]] = fz.flipMap()[j];
306 boolList subFlipStatus(subAddressing.size(), false);
307 forAll(subAddressing, j)
309 subFlipStatus[j] = fullFlipStatus[faceMap()[subAddressing[j]]];
312 fZonePtrs[i] = new faceZone
318 fvMeshSubsetPtr_->faceZones()
323 const cellZoneMesh& cellZones = baseMesh().cellZones();
325 List<cellZone*> cZonePtrs(cellZones.size());
329 const cellZone& cz = cellZones[i];
331 cZonePtrs[i] = new cellZone
334 subset(baseMesh().nCells(), cz, cellMap()),
336 fvMeshSubsetPtr_->cellZones()
342 fvMeshSubsetPtr_->addZones(pZonePtrs, fZonePtrs, cZonePtrs);
346 void Foam::fvMeshSubset::makeFvDictionaries()
348 // Try accessing dictionaries
353 this->name() + "Subset",
354 baseMesh().time().timeName(),
361 IOdictionary fvSchemes
373 if (!fvSchemes.headerOk())
375 Info << "Cannot read " << fvSchemes.path()
376 << ". Copy from base" << endl;
378 IOdictionary fvSchemesBase
383 baseMesh().time().system(),
390 fvSchemes = fvSchemesBase;
391 fvSchemes.regIOobject::write();
394 IOdictionary fvSolution
406 if (!fvSolution.headerOk())
408 Info << "Cannot read " << fvSolution.path()
409 << ". Copy from base" << endl;
411 IOdictionary fvSolutionBase
416 baseMesh().time().system(),
423 fvSolution = fvSolutionBase;
424 fvSolution.regIOobject::write();
429 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
431 // Construct from components
432 Foam::fvMeshSubset::fvMeshSubset
435 const fvMesh& baseMesh
440 fvMeshSubsetPtr_(NULL),
441 pointMeshSubsetPtr_(NULL),
449 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
451 Foam::fvMeshSubset::~fvMeshSubset()
453 deleteDemandDrivenData(fvMeshSubsetPtr_);
454 deleteDemandDrivenData(pointMeshSubsetPtr_);
458 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
460 void Foam::fvMeshSubset::setCellSubset
462 const labelHashSet& globalCellMap,
466 // Initial check on patches before doing anything time consuming.
467 const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
468 const cellList& oldCells = baseMesh().cells();
469 const faceList& oldFaces = baseMesh().allFaces();
470 const pointField& oldPoints = baseMesh().points();
471 const labelList& oldOwner = baseMesh().faceOwner();
472 const labelList& oldNeighbour = baseMesh().faceNeighbour();
474 label wantedPatchID = patchID;
476 if (wantedPatchID == -1)
478 // No explicit patch specified. Put in oldInternalFaces patch.
479 // Check if patch with this name already exists.
480 wantedPatchID = oldPatches.findPatchID("oldInternalFaces");
482 else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.size())
486 "fvMeshSubset::setCellSubset(const labelHashSet&"
487 ", const label patchID)"
488 ) << "Non-existing patch index " << wantedPatchID << endl
489 << "Should be between 0 and " << oldPatches.size()-1
490 << abort(FatalError);
494 cellMap_ = globalCellMap.toc();
496 // Sort the cell map in the ascending order
499 label nCellsInSet = cellMap_.size();
501 // Mark all used faces
503 Map<label> facesToSubset(primitiveMesh::facesPerCell_*nCellsInSet);
505 forAll (cellMap_, cellI)
507 // Mark all faces from the cell
508 const labelList& curFaces = oldCells[cellMap_[cellI]];
510 forAll (curFaces, faceI)
512 if (!facesToSubset.found(curFaces[faceI]))
514 facesToSubset.insert(curFaces[faceI], 1);
518 facesToSubset[curFaces[faceI]]++;
523 // Mark all used points and make a global-to-local face map
524 Map<label> globalFaceMap(facesToSubset.size());
526 // Make a global-to-local point map
527 Map<label> globalPointMap
529 primitiveMesh::pointsPerFace_*facesToSubset.size()
532 // This is done in two goes, so that the boundary faces are last
533 // in the list. Because of this, I need to create the face map
534 // along the way rather than just grab the table of contents.
535 labelList facesToc = facesToSubset.toc();
537 faceMap_.setSize(facesToc.size());
539 // 1. Get all faces that will be internal to the submesh.
540 forAll (facesToc, faceI)
542 if (facesToSubset[facesToc[faceI]] == 2)
544 // Mark face and increment number of points in set
545 faceMap_[globalFaceMap.size()] = facesToc[faceI];
546 globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
548 // Mark all points from the face
549 markPoints(oldFaces[facesToc[faceI]], globalPointMap);
553 // These are all the internal faces in the mesh.
554 label nInternalFaces = globalFaceMap.size();
557 // Where to insert old internal faces.
558 label oldPatchStart = labelMax;
559 if (wantedPatchID != -1)
561 oldPatchStart = oldPatches[wantedPatchID].start();
567 // 2. Boundary faces up to where we want to insert old internal faces
568 for (; faceI< facesToc.size(); faceI++)
570 if (facesToc[faceI] >= oldPatchStart)
576 !baseMesh().isInternalFace(facesToc[faceI])
577 && facesToSubset[facesToc[faceI]] == 1
580 // Mark face and increment number of points in set
581 faceMap_[globalFaceMap.size()] = facesToc[faceI];
582 globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
584 // Mark all points from the face
585 markPoints(oldFaces[facesToc[faceI]], globalPointMap);
589 // 3. old internal faces
590 forAll(facesToc, intFaceI)
594 baseMesh().isInternalFace(facesToc[intFaceI])
595 && facesToSubset[facesToc[intFaceI]] == 1
598 // Mark face and increment number of points in set
599 faceMap_[globalFaceMap.size()] = facesToc[intFaceI];
600 globalFaceMap.insert(facesToc[intFaceI], globalFaceMap.size());
602 // Mark all points from the face
603 markPoints(oldFaces[facesToc[intFaceI]], globalPointMap);
607 // 4. Remaining boundary faces
608 for (; faceI< facesToc.size(); faceI++)
612 !baseMesh().isInternalFace(facesToc[faceI])
613 && facesToSubset[facesToc[faceI]] == 1
616 // Mark face and increment number of points in set
617 faceMap_[globalFaceMap.size()] = facesToc[faceI];
618 globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
620 // Mark all points from the face
621 markPoints(oldFaces[facesToc[faceI]], globalPointMap);
627 // Grab the points map
628 pointMap_ = globalPointMap.toc();
631 forAll (pointMap_, pointI)
633 globalPointMap[pointMap_[pointI]] = pointI;
638 Pout<< "Number of cells in new mesh: " << nCellsInSet << nl
639 << "Number of faces in new mesh: " << globalFaceMap.size() << nl
640 << "Number of points in new mesh: " << globalPointMap.size()
645 pointField newPoints(globalPointMap.size());
647 label nNewPoints = 0;
649 forAll (pointMap_, pointI)
651 newPoints[nNewPoints] = oldPoints[pointMap_[pointI]];
655 faceList newFaces(globalFaceMap.size());
659 // Make internal faces
660 for (label faceI = 0; faceI < nInternalFaces; faceI++)
662 const face& oldF = oldFaces[faceMap_[faceI]];
664 face newF(oldF.size());
668 newF[i] = globalPointMap[oldF[i]];
671 newFaces[nNewFaces] = newF;
675 // Make boundary faces
677 label nbSize = oldPatches.size();
678 label oldInternalPatchID = -1;
680 if (wantedPatchID == -1)
682 // Create 'oldInternalFaces' patch at the end
683 // and put all exposed internal faces in there.
684 oldInternalPatchID = nbSize;
690 oldInternalPatchID = wantedPatchID;
694 // Grad size and start of each patch on the fly. Because of the
695 // structure of the underlying mesh, the patches will appear in the
697 labelList boundaryPatchSizes(nbSize, 0);
699 // Assign boundary faces. Visited in order of faceMap_.
700 for (label faceI = nInternalFaces; faceI < faceMap_.size(); faceI++)
702 label oldFaceI = faceMap_[faceI];
704 face oldF = oldFaces[oldFaceI];
706 // Turn the faces as necessary to point outwards
707 if (baseMesh().isInternalFace(oldFaceI))
709 // Internal face. Possibly turned the wrong way round
712 !globalCellMap.found(oldOwner[oldFaceI])
713 && globalCellMap.found(oldNeighbour[oldFaceI])
716 oldF = oldFaces[oldFaceI].reverseFace();
719 // Update count for patch
720 boundaryPatchSizes[oldInternalPatchID]++;
724 // Boundary face. Increment the appropriate patch
725 label patchOfFace = oldPatches.whichPatch(oldFaceI);
727 // Update count for patch
728 boundaryPatchSizes[patchOfFace]++;
731 face newF(oldF.size());
735 newF[i] = globalPointMap[oldF[i]];
738 newFaces[nNewFaces] = newF;
745 cellList newCells(nCellsInSet);
749 forAll (cellMap_, cellI)
751 const labelList& oldC = oldCells[cellMap_[cellI]];
753 labelList newC(oldC.size());
757 newC[i] = globalFaceMap[oldC[i]];
760 newCells[nNewCells] = cell(newC);
766 fvMeshSubsetPtr_ = new fvMesh
770 this->name() + "Subset",
771 baseMesh().time().timeName(),
782 deleteDemandDrivenData(pointMeshSubsetPtr_);
786 List<polyPatch*> newBoundary(nbSize);
787 patchMap_.setSize(nbSize);
788 label nNewPatches = 0;
789 label patchStart = nInternalFaces;
791 forAll(oldPatches, patchI)
793 if (boundaryPatchSizes[patchI] > 0)
795 // Patch still exists. Add it
796 newBoundary[nNewPatches] = oldPatches[patchI].clone
798 fvMeshSubsetPtr_->boundaryMesh(),
800 boundaryPatchSizes[patchI],
804 patchStart += boundaryPatchSizes[patchI];
805 patchMap_[nNewPatches] = patchI;
810 if (wantedPatchID == -1)
812 // Newly created patch so is at end. Check if any faces in it.
813 if (boundaryPatchSizes[oldInternalPatchID] > 0)
815 newBoundary[nNewPatches] = new polyPatch
818 boundaryPatchSizes[oldInternalPatchID],
821 fvMeshSubsetPtr_->boundaryMesh()
824 // The index for the first patch is -1 as it originates from
825 // the internal faces
826 patchMap_[nNewPatches] = -1;
831 // Reset the patch lists
832 newBoundary.setSize(nNewPatches);
833 patchMap_.setSize(nNewPatches);
836 fvMeshSubsetPtr_->addFvPatches(newBoundary);
838 // Subset and add any zones
843 void Foam::fvMeshSubset::setLargeCellSubset
845 const labelList& region,
846 const label currentRegion,
851 const cellList& oldCells = baseMesh().cells();
852 const faceList& oldFaces = baseMesh().faces();
853 const pointField& oldPoints = baseMesh().points();
854 const labelList& oldOwner = baseMesh().faceOwner();
855 const labelList& oldNeighbour = baseMesh().faceNeighbour();
856 const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
857 const label oldNInternalFaces = baseMesh().nInternalFaces();
861 if (region.size() != oldCells.size())
865 "fvMeshSubset::setCellSubset(const labelList&"
866 ", const label, const label, const bool)"
867 ) << "Size of region " << region.size()
868 << " is not equal to number of cells in mesh " << oldCells.size()
869 << abort(FatalError);
873 label wantedPatchID = patchID;
875 if (wantedPatchID == -1)
877 // No explicit patch specified. Put in oldInternalFaces patch.
878 // Check if patch with this name already exists.
879 wantedPatchID = oldPatches.findPatchID("oldInternalFaces");
881 else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.size())
885 "fvMeshSubset::setCellSubset(const labelList&"
886 ", const label, const label, const bool)"
887 ) << "Non-existing patch index " << wantedPatchID << endl
888 << "Should be between 0 and " << oldPatches.size()-1
889 << abort(FatalError);
893 // Get the cells for the current region.
894 cellMap_.setSize(oldCells.size());
895 label nCellsInSet = 0;
897 forAll (region, oldCellI)
899 if (region[oldCellI] == currentRegion)
901 cellMap_[nCellsInSet++] = oldCellI;
904 cellMap_.setSize(nCellsInSet);
907 // Mark all used faces. Count number of cells using them
908 // 0: face not used anymore
909 // 1: face used by one cell, face becomes/stays boundary face
910 // 2: face still used and remains internal face
911 // 3: face coupled and used by one cell only (so should become normal,
912 // non-coupled patch face)
914 // Note that this is not really nessecary - but means we can size things
915 // correctly. Also makes handling coupled faces much easier.
917 labelList nCellsUsingFace(oldFaces.size(), 0);
919 label nFacesInSet = 0;
920 forAll (oldFaces, oldFaceI)
922 bool faceUsed = false;
924 if (region[oldOwner[oldFaceI]] == currentRegion)
926 nCellsUsingFace[oldFaceI]++;
932 baseMesh().isInternalFace(oldFaceI)
933 && (region[oldNeighbour[oldFaceI]] == currentRegion)
936 nCellsUsingFace[oldFaceI]++;
945 faceMap_.setSize(nFacesInSet);
947 // Handle coupled faces. Modifies patch faces to be uncoupled to 3.
948 doCoupledPatches(syncPar, nCellsUsingFace);
951 // See which patch to use for exposed internal faces.
952 label oldInternalPatchID = 0;
954 // Insert faces before which patch
955 label nextPatchID = oldPatches.size();
957 // old to new patches
958 labelList globalPatchMap(oldPatches.size());
961 label nbSize = oldPatches.size();
963 if (wantedPatchID == -1)
965 // Create 'oldInternalFaces' patch at the end (or before
967 // and put all exposed internal faces in there.
969 forAll(oldPatches, patchI)
971 if (isA<processorPolyPatch>(oldPatches[patchI]))
973 nextPatchID = patchI;
976 oldInternalPatchID++;
981 // adapt old to new patches for inserted patch
982 for (label oldPatchI = 0; oldPatchI < nextPatchID; oldPatchI++)
984 globalPatchMap[oldPatchI] = oldPatchI;
988 label oldPatchI = nextPatchID;
989 oldPatchI < oldPatches.size();
993 globalPatchMap[oldPatchI] = oldPatchI+1;
998 oldInternalPatchID = wantedPatchID;
999 nextPatchID = wantedPatchID+1;
1001 // old to new patches
1002 globalPatchMap = identity(oldPatches.size());
1005 labelList boundaryPatchSizes(nbSize, 0);
1008 // Make a global-to-local point map
1009 labelList globalPointMap(oldPoints.size(), -1);
1011 labelList globalFaceMap(oldFaces.size(), -1);
1014 // 1. Pick up all preserved internal faces.
1015 for (label oldFaceI = 0; oldFaceI < oldNInternalFaces; oldFaceI++)
1017 if (nCellsUsingFace[oldFaceI] == 2)
1019 globalFaceMap[oldFaceI] = faceI;
1020 faceMap_[faceI++] = oldFaceI;
1022 // Mark all points from the face
1023 markPoints(oldFaces[oldFaceI], globalPointMap);
1027 // These are all the internal faces in the mesh.
1028 label nInternalFaces = faceI;
1030 // 2. Boundary faces up to where we want to insert old internal faces
1033 label oldPatchI = 0;
1034 oldPatchI < oldPatches.size()
1035 && oldPatchI < nextPatchID;
1039 const polyPatch& oldPatch = oldPatches[oldPatchI];
1041 label oldFaceI = oldPatch.start();
1043 forAll (oldPatch, i)
1045 if (nCellsUsingFace[oldFaceI] == 1)
1047 // Boundary face is kept.
1049 // Mark face and increment number of points in set
1050 globalFaceMap[oldFaceI] = faceI;
1051 faceMap_[faceI++] = oldFaceI;
1053 // Mark all points from the face
1054 markPoints(oldFaces[oldFaceI], globalPointMap);
1056 // Increment number of patch faces
1057 boundaryPatchSizes[globalPatchMap[oldPatchI]]++;
1063 // 3a. old internal faces that have become exposed.
1064 for (label oldFaceI = 0; oldFaceI < oldNInternalFaces; oldFaceI++)
1066 if (nCellsUsingFace[oldFaceI] == 1)
1068 globalFaceMap[oldFaceI] = faceI;
1069 faceMap_[faceI++] = oldFaceI;
1071 // Mark all points from the face
1072 markPoints(oldFaces[oldFaceI], globalPointMap);
1074 // Increment number of patch faces
1075 boundaryPatchSizes[oldInternalPatchID]++;
1079 // 3b. coupled patch faces that have become uncoupled.
1082 label oldFaceI = oldNInternalFaces;
1083 oldFaceI < oldFaces.size();
1087 if (nCellsUsingFace[oldFaceI] == 3)
1089 globalFaceMap[oldFaceI] = faceI;
1090 faceMap_[faceI++] = oldFaceI;
1092 // Mark all points from the face
1093 markPoints(oldFaces[oldFaceI], globalPointMap);
1095 // Increment number of patch faces
1096 boundaryPatchSizes[oldInternalPatchID]++;
1100 // 4. Remaining boundary faces
1103 label oldPatchI = nextPatchID;
1104 oldPatchI < oldPatches.size();
1108 const polyPatch& oldPatch = oldPatches[oldPatchI];
1110 label oldFaceI = oldPatch.start();
1112 forAll (oldPatch, i)
1114 if (nCellsUsingFace[oldFaceI] == 1)
1116 // Boundary face is kept.
1118 // Mark face and increment number of points in set
1119 globalFaceMap[oldFaceI] = faceI;
1120 faceMap_[faceI++] = oldFaceI;
1122 // Mark all points from the face
1123 markPoints(oldFaces[oldFaceI], globalPointMap);
1125 // Increment number of patch faces
1126 boundaryPatchSizes[globalPatchMap[oldPatchI]]++;
1132 if (faceI != nFacesInSet)
1136 "fvMeshSubset::setCellSubset(const labelList&"
1137 ", const label, const label, const bool)"
1138 ) << "Problem" << abort(FatalError);
1142 // Grab the points map
1143 label nPointsInSet = 0;
1145 forAll(globalPointMap, pointI)
1147 if (globalPointMap[pointI] != -1)
1152 pointMap_.setSize(nPointsInSet);
1156 forAll(globalPointMap, pointI)
1158 if (globalPointMap[pointI] != -1)
1160 pointMap_[nPointsInSet] = pointI;
1161 globalPointMap[pointI] = nPointsInSet;
1166 Pout<< "Number of cells in new mesh : " << cellMap_.size() << endl;
1167 Pout<< "Number of faces in new mesh : " << faceMap_.size() << endl;
1168 Pout<< "Number of points in new mesh: " << pointMap_.size() << endl;
1171 pointField newPoints(pointMap_.size());
1173 label nNewPoints = 0;
1175 forAll (pointMap_, pointI)
1177 newPoints[nNewPoints] = oldPoints[pointMap_[pointI]];
1181 faceList newFaces(faceMap_.size());
1183 label nNewFaces = 0;
1185 // Make internal faces
1186 for (label faceI = 0; faceI < nInternalFaces; faceI++)
1188 const face& oldF = oldFaces[faceMap_[faceI]];
1190 face newF(oldF.size());
1194 newF[i] = globalPointMap[oldF[i]];
1197 newFaces[nNewFaces] = newF;
1202 // Make boundary faces. (different from internal since might need to be
1204 for (label faceI = nInternalFaces; faceI < faceMap_.size(); faceI++)
1206 label oldFaceI = faceMap_[faceI];
1208 face oldF = oldFaces[oldFaceI];
1210 // Turn the faces as necessary to point outwards
1211 if (baseMesh().isInternalFace(oldFaceI))
1213 // Was internal face. Possibly turned the wrong way round
1216 region[oldOwner[oldFaceI]] != currentRegion
1217 && region[oldNeighbour[oldFaceI]] == currentRegion
1220 oldF = oldFaces[oldFaceI].reverseFace();
1224 // Relabel vertices of the (possibly turned) face.
1225 face newF(oldF.size());
1229 newF[i] = globalPointMap[oldF[i]];
1232 newFaces[nNewFaces] = newF;
1239 cellList newCells(nCellsInSet);
1241 label nNewCells = 0;
1243 forAll (cellMap_, cellI)
1245 const labelList& oldC = oldCells[cellMap_[cellI]];
1247 labelList newC(oldC.size());
1251 newC[i] = globalFaceMap[oldC[i]];
1254 newCells[nNewCells] = cell(newC);
1260 makeFvDictionaries();
1262 fvMeshSubsetPtr_ = new fvMesh
1266 this->name() + "Subset",
1267 baseMesh().time().timeName(),
1272 xferMove(newPoints),
1275 syncPar // parallel synchronisation
1279 deleteDemandDrivenData(pointMeshSubsetPtr_);
1282 List<polyPatch*> newBoundary(nbSize);
1283 patchMap_.setSize(nbSize);
1284 label nNewPatches = 0;
1285 label patchStart = nInternalFaces;
1288 // For parallel: only remove patch if none of the processors has it.
1289 // This only gets done for patches before the one being inserted
1290 // (so patches < nextPatchID)
1292 // Get sum of patch sizes. Zero if patch can be deleted.
1293 labelList globalPatchSizes(boundaryPatchSizes);
1294 globalPatchSizes.setSize(nextPatchID);
1296 if (syncPar && Pstream::parRun())
1298 // Get patch names (up to nextPatchID)
1299 List<wordList> patchNames(Pstream::nProcs());
1300 patchNames[Pstream::myProcNo()] = oldPatches.names();
1301 patchNames[Pstream::myProcNo()].setSize(nextPatchID);
1302 Pstream::gatherList(patchNames);
1303 Pstream::scatterList(patchNames);
1305 // Get patch sizes (up to nextPatchID).
1306 // Note that up to nextPatchID the globalPatchMap is an identity so
1307 // no need to index through that.
1308 Pstream::listCombineGather(globalPatchSizes, plusEqOp<label>());
1309 Pstream::listCombineScatter(globalPatchSizes);
1311 // Now all processors have all the patchnames.
1312 // Decide: if all processors have the same patch names and size is zero
1313 // everywhere remove the patch.
1314 bool samePatches = true;
1316 for (label procI = 1; procI < patchNames.size(); procI++)
1318 if (patchNames[procI] != patchNames[0])
1320 samePatches = false;
1327 // Patchnames not sync on all processors so disable removal of
1328 // zero sized patches.
1329 globalPatchSizes = labelMax;
1338 label oldPatchI = 0;
1339 oldPatchI < oldPatches.size()
1340 && oldPatchI < nextPatchID;
1344 label newSize = boundaryPatchSizes[globalPatchMap[oldPatchI]];
1346 // Clone (even if 0 size)
1347 newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
1349 fvMeshSubsetPtr_->boundaryMesh(),
1355 patchStart += newSize;
1356 patchMap_[nNewPatches] = oldPatchI; // compact patchMap
1362 if (wantedPatchID == -1)
1364 label oldInternalSize = boundaryPatchSizes[oldInternalPatchID];
1368 reduce(oldInternalSize, sumOp<label>());
1371 // Newly created patch so is at end. Check if any faces in it.
1372 if (oldInternalSize > 0)
1374 newBoundary[nNewPatches] = new polyPatch
1377 boundaryPatchSizes[oldInternalPatchID],
1380 fvMeshSubsetPtr_->boundaryMesh()
1383 Pout<< " oldInternalFaces : "
1384 << boundaryPatchSizes[oldInternalPatchID] << endl;
1386 // The index for the first patch is -1 as it originates from
1387 // the internal faces
1388 patchStart += boundaryPatchSizes[oldInternalPatchID];
1389 patchMap_[nNewPatches] = -1;
1398 label oldPatchI = nextPatchID;
1399 oldPatchI < oldPatches.size();
1403 label newSize = boundaryPatchSizes[globalPatchMap[oldPatchI]];
1405 // Patch still exists. Add it
1406 newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
1408 fvMeshSubsetPtr_->boundaryMesh(),
1414 //Pout<< " " << oldPatches[oldPatchI].name() << " : "
1415 // << newSize << endl;
1417 patchStart += newSize;
1418 patchMap_[nNewPatches] = oldPatchI; // compact patchMap
1423 // Reset the patch lists
1424 newBoundary.setSize(nNewPatches);
1425 patchMap_.setSize(nNewPatches);
1428 // Add the fvPatches
1429 fvMeshSubsetPtr_->addFvPatches(newBoundary, syncPar);
1431 // Subset and add any zones
1436 void Foam::fvMeshSubset::setLargeCellSubset
1438 const labelHashSet& globalCellMap,
1439 const label patchID,
1443 labelList region(baseMesh().nCells(), 0);
1445 forAllConstIter (labelHashSet, globalCellMap, iter)
1447 region[iter.key()] = 1;
1449 setLargeCellSubset(region, 1, patchID, syncPar);
1453 const fvMesh& Foam::fvMeshSubset::subMesh() const
1457 return *fvMeshSubsetPtr_;
1461 fvMesh& Foam::fvMeshSubset::subMesh()
1465 return *fvMeshSubsetPtr_;
1469 const pointMesh& Foam::fvMeshSubset::subPointMesh() const
1471 if (!pointMeshSubsetPtr_)
1473 pointMeshSubsetPtr_ = new pointMesh(subMesh());
1476 return *pointMeshSubsetPtr_;
1480 pointMesh& Foam::fvMeshSubset::subPointMesh()
1482 if (!pointMeshSubsetPtr_)
1484 pointMeshSubsetPtr_ = new pointMesh(subMesh());
1487 return *pointMeshSubsetPtr_;
1491 const labelList& Foam::fvMeshSubset::pointMap() const
1499 const labelList& Foam::fvMeshSubset::faceMap() const
1507 const labelList& Foam::fvMeshSubset::cellMap() const
1515 const labelList& Foam::fvMeshSubset::patchMap() const
1523 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1525 } // End namespace Foam
1527 // ************************************************************************* //