1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
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
25 \*---------------------------------------------------------------------------*/
27 #include "polyTopoChange.H"
28 #include "SortableList.H"
30 #include "polyAddPoint.H"
31 #include "polyModifyPoint.H"
32 #include "polyRemovePoint.H"
33 #include "polyAddFace.H"
34 #include "polyModifyFace.H"
35 #include "polyRemoveFace.H"
36 #include "polyAddCell.H"
37 #include "polyModifyCell.H"
38 #include "polyRemoveCell.H"
39 #include "objectMap.H"
40 #include "processorPolyPatch.H"
43 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
47 defineTypeNameAndDebug(polyTopoChange, 0);
51 const Foam::point Foam::polyTopoChange::greatPoint
59 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
61 // Renumber with special handling for merged items (marked with <-1)
62 void Foam::polyTopoChange::renumberReverseMap
65 DynamicList<label>& elems
70 label val = elems[elemI];
74 elems[elemI] = map[val];
78 label mergedVal = -val-2;
79 elems[elemI] = -map[mergedVal]-2;
85 void Foam::polyTopoChange::renumber
91 labelHashSet newElems(elems.size());
93 forAllConstIter(labelHashSet, elems, iter)
95 label newElem = map[iter.key()];
99 newElems.insert(newElem);
103 elems.transfer(newElems);
107 // Renumber and remove -1 elements.
108 void Foam::polyTopoChange::renumberCompact
110 const labelList& map,
118 label newVal = map[elems[elemI]];
122 elems[newElemI++] = newVal;
125 elems.setSize(newElemI);
129 void Foam::polyTopoChange::countMap
131 const labelList& map,
132 const labelList& reverseMap,
144 forAll(map, newCellI)
146 label oldCellI = map[newCellI];
150 if (reverseMap[oldCellI] == newCellI)
156 // Added (from another cell v.s. inflated from face/point)
160 else if (oldCellI == -1)
162 // Created from nothing
167 FatalErrorIn("countMap") << "old:" << oldCellI
168 << " new:" << newCellI << abort(FatalError);
172 forAll(reverseMap, oldCellI)
174 label newCellI = reverseMap[oldCellI];
180 else if (newCellI == -1)
187 // merged into -newCellI-2
194 Foam::labelHashSet Foam::polyTopoChange::getSetIndices
196 const PackedBoolList& lst
199 labelHashSet values(lst.count());
211 void Foam::polyTopoChange::writeMeshStats(const polyMesh& mesh, Ostream& os)
213 const polyBoundaryMesh& patches = mesh.boundaryMesh();
215 labelList patchSizes(patches.size());
216 labelList patchStarts(patches.size());
217 forAll(patches, patchI)
219 patchSizes[patchI] = patches[patchI].size();
220 patchStarts[patchI] = patches[patchI].start();
223 os << " Points : " << mesh.nPoints() << nl
224 << " Faces : " << mesh.nFaces() << nl
225 << " Cells : " << mesh.nCells() << nl
226 << " PatchSizes : " << patchSizes << nl
227 << " PatchStarts : " << patchStarts << nl
232 void Foam::polyTopoChange::getMergeSets
234 const labelList& reverseCellMap,
235 const labelList& cellMap,
236 List<objectMap>& cellsFromCells
239 // Per new cell the number of old cells that have been merged into it
240 labelList nMerged(cellMap.size(), 1);
242 forAll(reverseCellMap, oldCellI)
244 label newCellI = reverseCellMap[oldCellI];
248 label mergeCellI = -newCellI-2;
250 nMerged[mergeCellI]++;
254 // From merged cell to set index
255 labelList cellToMergeSet(cellMap.size(), -1);
259 forAll(nMerged, cellI)
261 if (nMerged[cellI] > 1)
263 cellToMergeSet[cellI] = nSets++;
267 // Collect cell labels.
268 // Each objectMap will have
269 // - index : new mesh cell label
270 // - masterObjects : list of old cells that have been merged. Element 0
271 // will be the original destination cell label.
273 cellsFromCells.setSize(nSets);
275 forAll(reverseCellMap, oldCellI)
277 label newCellI = reverseCellMap[oldCellI];
281 label mergeCellI = -newCellI-2;
283 // oldCellI was merged into mergeCellI
285 label setI = cellToMergeSet[mergeCellI];
287 objectMap& mergeSet = cellsFromCells[setI];
289 if (mergeSet.masterObjects().empty())
291 // First occurrence of master cell mergeCellI
293 mergeSet.index() = mergeCellI;
294 mergeSet.masterObjects().setSize(nMerged[mergeCellI]);
297 mergeSet.masterObjects()[0] = cellMap[mergeCellI];
300 mergeSet.masterObjects()[1] = oldCellI;
302 nMerged[mergeCellI] = 2;
306 mergeSet.masterObjects()[nMerged[mergeCellI]++] = oldCellI;
313 void Foam::polyTopoChange::checkFace
325 if (own == -1 && zoneI != -1)
329 else if (patchI == -1 || patchI >= nPatches_)
333 "polyTopoChange::checkFace(const face&, const label"
334 ", const label, const label, const label)"
335 ) << "Face has no neighbour (so external) but does not have"
336 << " a valid patch" << nl
338 << " faceI(-1 if added face):" << faceI
339 << " own:" << own << " nei:" << nei
340 << " patchI:" << patchI << abort(FatalError);
349 "polyTopoChange::checkFace(const face&, const label"
350 ", const label, const label, const label)"
351 ) << "Cannot both have valid patchI and neighbour" << nl
353 << " faceI(-1 if added face):" << faceI
354 << " own:" << own << " nei:" << nei
355 << " patchI:" << patchI << abort(FatalError);
362 "polyTopoChange::checkFace(const face&, const label"
363 ", const label, const label, const label)"
364 ) << "Owner cell label should be less than neighbour cell label"
367 << " faceI(-1 if added face):" << faceI
368 << " own:" << own << " nei:" << nei
369 << " patchI:" << patchI << abort(FatalError);
373 if (f.size() < 3 || findIndex(f, -1) != -1)
377 "polyTopoChange::checkFace(const face&, const label"
378 ", const label, const label, const label)"
379 ) << "Illegal vertices in face"
382 << " faceI(-1 if added face):" << faceI
383 << " own:" << own << " nei:" << nei
384 << " patchI:" << patchI << abort(FatalError);
386 if (faceI >= 0 && faceI < faces_.size() && faceRemoved(faceI))
390 "polyTopoChange::checkFace(const face&, const label"
391 ", const label, const label, const label)"
392 ) << "Face already marked for removal"
395 << " faceI(-1 if added face):" << faceI
396 << " own:" << own << " nei:" << nei
397 << " patchI:" << patchI << abort(FatalError);
401 if (f[fp] < points_.size() && pointRemoved(f[fp]))
405 "polyTopoChange::checkFace(const face&, const label"
406 ", const label, const label, const label)"
407 ) << "Face uses removed vertices"
410 << " faceI(-1 if added face):" << faceI
411 << " own:" << own << " nei:" << nei
412 << " patchI:" << patchI << abort(FatalError);
418 void Foam::polyTopoChange::makeCells
420 const label nActiveFaces,
421 labelList& cellFaces,
422 labelList& cellFaceOffsets
425 cellFaces.setSize(2*nActiveFaces);
426 cellFaceOffsets.setSize(cellMap_.size() + 1);
429 labelList nNbrs(cellMap_.size(), 0);
431 // 1. Count faces per cell
433 for (label faceI = 0; faceI < nActiveFaces; faceI++)
435 nNbrs[faceOwner_[faceI]]++;
437 for (label faceI = 0; faceI < nActiveFaces; faceI++)
439 if (faceNeighbour_[faceI] >= 0)
441 nNbrs[faceNeighbour_[faceI]]++;
445 // 2. Calculate offsets
447 cellFaceOffsets[0] = 0;
450 cellFaceOffsets[cellI+1] = cellFaceOffsets[cellI] + nNbrs[cellI];
453 // 3. Fill faces per cell
455 // reset the whole list to use as counter
458 for (label faceI = 0; faceI < nActiveFaces; faceI++)
460 label cellI = faceOwner_[faceI];
462 cellFaces[cellFaceOffsets[cellI] + nNbrs[cellI]++] = faceI;
465 for (label faceI = 0; faceI < nActiveFaces; faceI++)
467 label cellI = faceNeighbour_[faceI];
471 cellFaces[cellFaceOffsets[cellI] + nNbrs[cellI]++] = faceI;
475 // Last offset points to beyond end of cellFaces.
476 cellFaces.setSize(cellFaceOffsets[cellMap_.size()]);
480 // Create cell-cell addressing. Called after compaction (but before ordering)
482 void Foam::polyTopoChange::makeCellCells
484 const label nActiveFaces,
485 CompactListList<label>& cellCells
488 // Neighbours per cell
489 labelList nNbrs(cellMap_.size(), 0);
491 // Overall number of cellCells
492 label nCellCells = 0;
494 // 1. Count neighbours (through internal faces) per cell
496 for (label faceI = 0; faceI < nActiveFaces; faceI++)
498 if (faceNeighbour_[faceI] >= 0)
500 nNbrs[faceOwner_[faceI]]++;
501 nNbrs[faceNeighbour_[faceI]]++;
506 cellCells.setSize(cellMap_.size(), nCellCells);
508 // 2. Calculate offsets
510 labelList& offsets = cellCells.offsets();
515 sumSize += nNbrs[cellI];
516 offsets[cellI] = sumSize;
519 // 3. Fill faces per cell
521 // reset the whole list to use as counter
524 for (label faceI = 0; faceI < nActiveFaces; faceI++)
526 label nei = faceNeighbour_[faceI];
530 label own = faceOwner_[faceI];
531 cellCells.m()[cellCells.index(own, nNbrs[own]++)] = nei;
532 cellCells.m()[cellCells.index(nei, nNbrs[nei]++)] = own;
538 // Cell ordering (based on bandCompression).
539 // Handles removed cells. Returns number of remaining cells.
540 Foam::label Foam::polyTopoChange::getCellOrder
542 const CompactListList<label>& cellCellAddressing,
546 const labelList& offsets = cellCellAddressing.offsets();
548 labelList newOrder(cellCellAddressing.size());
550 // Fifo buffer for string of cells
551 SLList<label> nextCell;
553 // Whether cell has been done already
554 PackedBoolList visited(cellCellAddressing.size());
556 label cellInOrder = 0;
559 // loop over the cells
560 forAll (visited, cellI)
562 // find the first non-removed cell that has not been visited yet
563 if (!cellRemoved(cellI) && visited.get(cellI) == 0)
565 // use this cell as a start
566 nextCell.append(cellI);
568 // loop through the nextCell list. Add the first cell into the
569 // cell order if it has not already been visited and ask for its
570 // neighbours. If the neighbour in question has not been visited,
571 // add it to the end of the nextCell list
575 label currentCell = nextCell.removeHead();
577 if (visited.get(currentCell) == 0)
579 visited.set(currentCell, 1);
581 // add into cellOrder
582 newOrder[cellInOrder] = currentCell;
585 // find if the neighbours have been visited
586 label i0 = (currentCell == 0 ? 0 : offsets[currentCell-1]);
587 label i1 = offsets[currentCell];
589 for (label i = i0; i < i1; i++)
591 label nbr = cellCellAddressing.m()[i];
593 if (!cellRemoved(nbr) && visited.get(nbr) == 0)
595 // not visited, add to the list
596 nextCell.append(nbr);
601 while (nextCell.size());
606 // Now we have new-to-old in newOrder.
607 newOrder.setSize(cellInOrder);
609 // Invert to get old-to-new. Make sure removed (i.e. unmapped) cells are -1.
610 oldToNew = invert(cellCellAddressing.size(), newOrder);
616 // Determine order for faces:
617 // - upper-triangular order for internal faces
618 // - external faces after internal faces and in patch order.
619 void Foam::polyTopoChange::getFaceOrder
621 const label nActiveFaces,
622 const labelList& cellFaces,
623 const labelList& cellFaceOffsets,
626 labelList& patchSizes,
627 labelList& patchStarts
630 oldToNew.setSize(faceOwner_.size());
633 // First unassigned face
636 forAll(cellMap_, cellI)
638 label startOfCell = cellFaceOffsets[cellI];
639 label nFaces = cellFaceOffsets[cellI+1] - startOfCell;
641 // Neighbouring cells
642 SortableList<label> nbr(nFaces);
644 for (label i = 0; i < nFaces; i++)
646 label faceI = cellFaces[startOfCell + i];
648 label nbrCellI = faceNeighbour_[faceI];
650 if (faceI >= nActiveFaces)
655 else if (nbrCellI != -1)
657 // Internal face. Get cell on other side.
658 if (nbrCellI == cellI)
660 nbrCellI = faceOwner_[faceI];
663 if (cellI < nbrCellI)
670 // nbrCell is master. Let it handle this face.
676 // External face. Do later.
687 oldToNew[cellFaces[startOfCell + nbr.indices()[i]]] =
694 // Pick up all patch faces in patch face order.
695 patchStarts.setSize(nPatches_);
697 patchSizes.setSize(nPatches_);
700 patchStarts[0] = newFaceI;
702 for (label faceI = 0; faceI < nActiveFaces; faceI++)
704 if (region_[faceI] >= 0)
706 patchSizes[region_[faceI]]++;
710 label faceI = patchStarts[0];
712 forAll(patchStarts, patchI)
714 patchStarts[patchI] = faceI;
715 faceI += patchSizes[patchI];
720 // Pout<< "patchSizes:" << patchSizes << nl
721 // << "patchStarts:" << patchStarts << endl;
724 labelList workPatchStarts(patchStarts);
726 for (label faceI = 0; faceI < nActiveFaces; faceI++)
728 if (region_[faceI] >= 0)
730 oldToNew[faceI] = workPatchStarts[region_[faceI]]++;
735 for (label faceI = nActiveFaces; faceI < oldToNew.size(); faceI++)
737 oldToNew[faceI] = faceI;
740 // Check done all faces.
741 forAll(oldToNew, faceI)
743 if (oldToNew[faceI] == -1)
747 "polyTopoChange::getFaceOrder"
748 "(const label, const labelList&, const labelList&)"
750 ) << "Did not determine new position"
751 << " for face " << faceI
752 << abort(FatalError);
758 // Reorder and compact faces according to map.
759 void Foam::polyTopoChange::reorderCompactFaces
762 const labelList& oldToNew
765 reorder(oldToNew, faces_);
766 faces_.setCapacity(newSize);
768 reorder(oldToNew, region_);
769 region_.setCapacity(newSize);
771 reorder(oldToNew, faceOwner_);
772 faceOwner_.setCapacity(newSize);
774 reorder(oldToNew, faceNeighbour_);
775 faceNeighbour_.setCapacity(newSize);
778 reorder(oldToNew, faceMap_);
779 faceMap_.setCapacity(newSize);
781 renumberReverseMap(oldToNew, reverseFaceMap_);
783 renumberKey(oldToNew, faceFromPoint_);
784 renumberKey(oldToNew, faceFromEdge_);
785 inplaceReorder(oldToNew, flipFaceFlux_);
786 flipFaceFlux_.setCapacity(newSize);
788 reorder(oldToNew, faceZone_);
789 faceZone_.setCapacity(newSize);
791 inplaceReorder(oldToNew, faceZoneFlip_);
792 faceZoneFlip_.setCapacity(newSize);
796 // Compact all and orders points and faces:
797 // - points into internal followed by external points
798 // - internalfaces upper-triangular
799 // - externalfaces after internal ones.
800 void Foam::polyTopoChange::compact
802 const bool orderCells,
803 const bool orderPoints,
804 label& nInternalPoints,
805 labelList& patchSizes,
806 labelList& patchStarts
811 reversePointMap_.shrink();
817 faceNeighbour_.shrink();
819 reverseFaceMap_.shrink();
823 reverseCellMap_.shrink();
828 label nActivePoints = 0;
830 labelList localPointMap(points_.size(), -1);
835 nInternalPoints = -1;
837 forAll(points_, pointI)
839 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
841 localPointMap[pointI] = newPointI++;
844 nActivePoints = newPointI;
848 forAll(points_, pointI)
850 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
856 // Mark boundary points
857 forAll(faceOwner_, faceI)
862 && faceOwner_[faceI] >= 0
863 && faceNeighbour_[faceI] < 0
866 // Valid boundary face
867 const face& f = faces_[faceI];
871 label pointI = f[fp];
873 if (localPointMap[pointI] == -1)
878 || retiredPoints_.found(pointI)
881 FatalErrorIn("polyTopoChange::compact(..)")
882 << "Removed or retired point " << pointI
884 << " at position " << faceI << endl
885 << "Probably face has not been adapted for"
886 << " removed points." << abort(FatalError);
888 localPointMap[pointI] = newPointI++;
894 label nBoundaryPoints = newPointI;
895 nInternalPoints = nActivePoints - nBoundaryPoints;
897 // Move the boundary addressing up
898 forAll(localPointMap, pointI)
900 if (localPointMap[pointI] != -1)
902 localPointMap[pointI] += nInternalPoints;
908 // Mark internal points
909 forAll(faceOwner_, faceI)
914 && faceOwner_[faceI] >= 0
915 && faceNeighbour_[faceI] >= 0
918 // Valid internal face
919 const face& f = faces_[faceI];
923 label pointI = f[fp];
925 if (localPointMap[pointI] == -1)
930 || retiredPoints_.found(pointI)
933 FatalErrorIn("polyTopoChange::compact(..)")
934 << "Removed or retired point " << pointI
936 << " at position " << faceI << endl
937 << "Probably face has not been adapted for"
938 << " removed points." << abort(FatalError);
940 localPointMap[pointI] = newPointI++;
946 if (newPointI != nInternalPoints)
948 FatalErrorIn("polyTopoChange::compact(..)")
949 << "Problem." << abort(FatalError);
951 newPointI = nActivePoints;
954 forAllConstIter(labelHashSet, retiredPoints_, iter)
956 localPointMap[iter.key()] = newPointI++;
962 Pout<< "Points : active:" << nActivePoints
963 << " removed:" << points_.size()-newPointI << endl;
966 reorder(localPointMap, points_);
967 points_.setCapacity(newPointI);
970 reorder(localPointMap, pointMap_);
971 pointMap_.setCapacity(newPointI);
972 renumberReverseMap(localPointMap, reversePointMap_);
974 reorder(localPointMap, pointZone_);
975 pointZone_.setCapacity(newPointI);
976 renumber(localPointMap, retiredPoints_);
978 // Use map to relabel face vertices
979 forAll(faces_, faceI)
981 face& f = faces_[faceI];
984 renumberCompact(localPointMap, f);
986 if (!faceRemoved(faceI) && f.size() < 3)
988 FatalErrorIn("polyTopoChange::compact(..)")
989 << "Created illegal face " << f
990 //<< " from face " << oldF
991 << " at position:" << faceI
992 << " when filtering removed points"
993 << abort(FatalError);
1001 labelList localFaceMap(faces_.size(), -1);
1004 forAll(faces_, faceI)
1006 if (!faceRemoved(faceI) && faceOwner_[faceI] >= 0)
1008 localFaceMap[faceI] = newFaceI++;
1011 nActiveFaces_ = newFaceI;
1013 forAll(faces_, faceI)
1015 if (!faceRemoved(faceI) && faceOwner_[faceI] < 0)
1018 localFaceMap[faceI] = newFaceI++;
1024 Pout<< "Faces : active:" << nActiveFaces_
1025 << " removed:" << faces_.size()-newFaceI << endl;
1029 reorderCompactFaces(newFaceI, localFaceMap);
1034 labelList localCellMap;
1039 // Construct cellCell addressing
1040 CompactListList<label> cellCells;
1041 makeCellCells(nActiveFaces_, cellCells);
1043 // Cell ordering (based on bandCompression). Handles removed cells.
1044 newCellI = getCellOrder(cellCells, localCellMap);
1048 // Compact out removed cells
1049 localCellMap.setSize(cellMap_.size());
1053 forAll(cellMap_, cellI)
1055 if (!cellRemoved(cellI))
1057 localCellMap[cellI] = newCellI++;
1064 Pout<< "Cells : active:" << newCellI
1065 << " removed:" << cellMap_.size()-newCellI << endl;
1068 // Renumber -if cells reordered or -if cells removed
1069 if (orderCells || (newCellI != cellMap_.size()))
1071 reorder(localCellMap, cellMap_);
1072 cellMap_.setCapacity(newCellI);
1073 renumberReverseMap(localCellMap, reverseCellMap_);
1075 reorder(localCellMap, cellZone_);
1076 cellZone_.setCapacity(newCellI);
1078 renumberKey(localCellMap, cellFromPoint_);
1079 renumberKey(localCellMap, cellFromEdge_);
1080 renumberKey(localCellMap, cellFromFace_);
1082 // Renumber owner/neighbour. Take into account if neighbour suddenly
1083 // gets lower cell than owner.
1084 forAll(faceOwner_, faceI)
1086 label own = faceOwner_[faceI];
1087 label nei = faceNeighbour_[faceI];
1092 faceOwner_[faceI] = localCellMap[own];
1096 // Update neighbour.
1097 faceNeighbour_[faceI] = localCellMap[nei];
1099 // Check if face needs reversing.
1102 faceNeighbour_[faceI] >= 0
1103 && faceNeighbour_[faceI] < faceOwner_[faceI]
1106 faces_[faceI] = faces_[faceI].reverseFace();
1107 Swap(faceOwner_[faceI], faceNeighbour_[faceI]);
1108 flipFaceFlux_[faceI] =
1110 flipFaceFlux_[faceI]
1114 faceZoneFlip_[faceI] =
1116 faceZoneFlip_[faceI]
1125 // Update neighbour.
1126 faceNeighbour_[faceI] = localCellMap[nei];
1132 // Reorder faces into upper-triangular and patch ordering
1134 // Create cells (packed storage)
1135 labelList cellFaces;
1136 labelList cellFaceOffsets;
1137 makeCells(nActiveFaces_, cellFaces, cellFaceOffsets);
1139 // Do upper triangular order and patch sorting
1140 labelList localFaceMap;
1153 reorderCompactFaces(localFaceMap.size(), localFaceMap);
1158 // Find faces to interpolate to create value for new face. Only used if
1159 // face was inflated from edge or point. Internal faces should only be
1160 // created from internal faces, external faces only from external faces
1161 // (and ideally the same patch)
1162 // Is bit problematic if there are no faces to select, i.e. in polyDualMesh
1163 // an internal face can be created from a boundary edge with no internal
1164 // faces connected to it.
1165 Foam::labelList Foam::polyTopoChange::selectFaces
1167 const primitiveMesh& mesh,
1168 const labelList& faceLabels,
1169 const bool internalFacesOnly
1174 forAll(faceLabels, i)
1176 label faceI = faceLabels[i];
1178 if (internalFacesOnly == mesh.isInternalFace(faceI))
1184 labelList collectedFaces;
1188 // Did not find any faces of the correct type so just use any old
1190 collectedFaces = faceLabels;
1194 collectedFaces.setSize(nFaces);
1198 forAll(faceLabels, i)
1200 label faceI = faceLabels[i];
1202 if (internalFacesOnly == mesh.isInternalFace(faceI))
1204 collectedFaces[nFaces++] = faceI;
1209 return collectedFaces;
1213 // Calculate pointMap per patch (so from patch point label to old patch point
1215 void Foam::polyTopoChange::calcPatchPointMap
1217 const List<Map<label> >& oldPatchMeshPointMaps,
1218 const polyBoundaryMesh& boundary,
1219 labelListList& patchPointMap
1222 patchPointMap.setSize(boundary.size());
1224 forAll(boundary, patchI)
1226 const labelList& meshPoints = boundary[patchI].meshPoints();
1228 const Map<label>& oldMeshPointMap = oldPatchMeshPointMaps[patchI];
1230 labelList& curPatchPointRnb = patchPointMap[patchI];
1232 curPatchPointRnb.setSize(meshPoints.size());
1234 forAll(meshPoints, i)
1236 if (meshPoints[i] < pointMap_.size())
1238 // Check if old point was part of same patch
1239 Map<label>::const_iterator ozmpmIter = oldMeshPointMap.find
1241 pointMap_[meshPoints[i]]
1244 if (ozmpmIter != oldMeshPointMap.end())
1246 curPatchPointRnb[i] = ozmpmIter();
1250 curPatchPointRnb[i] = -1;
1255 curPatchPointRnb[i] = -1;
1262 void Foam::polyTopoChange::calcFaceInflationMaps
1264 const polyMesh& mesh,
1265 List<objectMap>& facesFromPoints,
1266 List<objectMap>& facesFromEdges,
1267 List<objectMap>& facesFromFaces
1270 // Faces inflated from points
1271 // ~~~~~~~~~~~~~~~~~~~~~~~~~~
1273 facesFromPoints.setSize(faceFromPoint_.size());
1275 if (faceFromPoint_.size())
1277 label nFacesFromPoints = 0;
1279 // Collect all still existing faces connected to this point.
1280 forAllConstIter(Map<label>, faceFromPoint_, iter)
1282 label newFaceI = iter.key();
1284 if (region_[newFaceI] == -1)
1286 // Get internal faces using point on old mesh
1287 facesFromPoints[nFacesFromPoints++] = objectMap
1293 mesh.pointFaces()[iter()],
1300 // Get patch faces using point on old mesh
1301 facesFromPoints[nFacesFromPoints++] = objectMap
1307 mesh.pointFaces()[iter()],
1316 // Faces inflated from edges
1317 // ~~~~~~~~~~~~~~~~~~~~~~~~~
1319 facesFromEdges.setSize(faceFromEdge_.size());
1321 if (faceFromEdge_.size())
1323 label nFacesFromEdges = 0;
1325 // Collect all still existing faces connected to this edge.
1326 forAllConstIter(Map<label>, faceFromEdge_, iter)
1328 label newFaceI = iter.key();
1330 if (region_[newFaceI] == -1)
1332 // Get internal faces using edge on old mesh
1333 facesFromEdges[nFacesFromEdges++] = objectMap
1339 mesh.edgeFaces(iter()),
1346 // Get patch faces using edge on old mesh
1347 facesFromEdges[nFacesFromEdges++] = objectMap
1353 mesh.edgeFaces(iter()),
1362 // Faces from face merging
1363 // ~~~~~~~~~~~~~~~~~~~~~~~
1374 void Foam::polyTopoChange::calcCellInflationMaps
1376 const polyMesh& mesh,
1377 List<objectMap>& cellsFromPoints,
1378 List<objectMap>& cellsFromEdges,
1379 List<objectMap>& cellsFromFaces,
1380 List<objectMap>& cellsFromCells
1383 cellsFromPoints.setSize(cellFromPoint_.size());
1385 if (cellFromPoint_.size())
1387 label nCellsFromPoints = 0;
1389 // Collect all still existing faces connected to this point.
1390 forAllConstIter(Map<label>, cellFromPoint_, iter)
1392 cellsFromPoints[nCellsFromPoints++] = objectMap
1395 mesh.pointCells()[iter()]
1401 cellsFromEdges.setSize(cellFromEdge_.size());
1403 if (cellFromEdge_.size())
1405 label nCellsFromEdges = 0;
1407 // Collect all still existing faces connected to this point.
1408 forAllConstIter(Map<label>, cellFromEdge_, iter)
1410 cellsFromEdges[nCellsFromEdges++] = objectMap
1413 mesh.edgeCells()[iter()]
1419 cellsFromFaces.setSize(cellFromFace_.size());
1421 if (cellFromFace_.size())
1423 label nCellsFromFaces = 0;
1425 labelList cellsAroundFace(2, -1);
1427 // Collect all still existing faces connected to this point.
1428 forAllConstIter(Map<label>, cellFromFace_, iter)
1430 label oldFaceI = iter();
1432 if (mesh.isInternalFace(oldFaceI))
1434 cellsAroundFace[0] = mesh.faceOwner()[oldFaceI];
1435 cellsAroundFace[1] = mesh.faceNeighbour()[oldFaceI];
1439 cellsAroundFace[0] = mesh.faceOwner()[oldFaceI];
1440 cellsAroundFace[1] = -1;
1443 cellsFromFaces[nCellsFromFaces++] = objectMap
1452 // Cells from cell merging
1453 // ~~~~~~~~~~~~~~~~~~~~~~~
1464 void Foam::polyTopoChange::resetZones
1466 const polyMesh& mesh,
1468 labelListList& pointZoneMap,
1469 labelListList& faceZoneFaceMap,
1470 labelListList& cellZoneMap
1476 pointZoneMap.setSize(mesh.pointZones().size());
1478 const pointZoneMesh& pointZones = mesh.pointZones();
1480 // Count points per zone
1482 labelList nPoints(pointZones.size(), 0);
1484 forAll(pointZone_, pointI)
1486 label zoneI = pointZone_[pointI];
1488 if (zoneI >= pointZones.size())
1492 "resetZones(const polyMesh&, polyMesh&, labelListList&"
1493 "labelListList&, labelListList&)"
1494 ) << "Illegal zoneID " << zoneI << " for point "
1495 << pointI << " coord " << mesh.points()[pointI]
1496 << abort(FatalError);
1505 // Distribute points per zone
1507 labelListList addressing(pointZones.size());
1508 forAll(addressing, zoneI)
1510 addressing[zoneI].setSize(nPoints[zoneI]);
1514 forAll(pointZone_, pointI)
1516 label zoneI = pointZone_[pointI];
1519 addressing[zoneI][nPoints[zoneI]++] = pointI;
1522 // Sort the addressing
1523 forAll(addressing, zoneI)
1525 stableSort(addressing[zoneI]);
1528 // So now we both have old zones and the new addressing.
1529 // Invert the addressing to get pointZoneMap.
1530 forAll(addressing, zoneI)
1532 const pointZone& oldZone = pointZones[zoneI];
1533 const labelList& newZoneAddr = addressing[zoneI];
1535 labelList& curPzRnb = pointZoneMap[zoneI];
1536 curPzRnb.setSize(newZoneAddr.size());
1538 forAll(newZoneAddr, i)
1540 if (newZoneAddr[i] < pointMap_.size())
1542 curPzRnb[i] = oldZone.whichPoint(pointMap_[newZoneAddr[i]]);
1551 // Reset the addresing on the zone
1552 newMesh.pointZones().clearAddressing();
1553 forAll(newMesh.pointZones(), zoneI)
1557 Pout<< "pointZone:" << zoneI
1558 << " name:" << newMesh.pointZones()[zoneI].name()
1559 << " size:" << addressing[zoneI].size()
1563 newMesh.pointZones()[zoneI] = addressing[zoneI];
1571 faceZoneFaceMap.setSize(mesh.faceZones().size());
1573 const faceZoneMesh& faceZones = mesh.faceZones();
1575 labelList nFaces(faceZones.size(), 0);
1577 forAll(faceZone_, faceI)
1579 label zoneI = faceZone_[faceI];
1581 if (zoneI >= faceZones.size())
1585 "resetZones(const polyMesh&, polyMesh&, labelListList&"
1586 "labelListList&, labelListList&)"
1587 ) << "Illegal zoneID " << zoneI << " for face "
1589 << abort(FatalError);
1597 labelListList addressing(faceZones.size());
1598 boolListList flipMode(faceZones.size());
1600 forAll(addressing, zoneI)
1602 addressing[zoneI].setSize(nFaces[zoneI]);
1603 flipMode[zoneI].setSize(nFaces[zoneI]);
1607 forAll(faceZone_, faceI)
1609 label zoneI = faceZone_[faceI];
1613 label index = nFaces[zoneI]++;
1614 addressing[zoneI][index] = faceI;
1615 flipMode[zoneI][index] = faceZoneFlip_[faceI];
1618 // Sort the addressing
1619 forAll(addressing, zoneI)
1622 sortedOrder(addressing[zoneI], newToOld);
1624 labelList newAddressing(addressing[zoneI].size());
1625 forAll(newAddressing, i)
1627 newAddressing[i] = addressing[zoneI][newToOld[i]];
1629 addressing[zoneI].transfer(newAddressing);
1632 boolList newFlipMode(flipMode[zoneI].size());
1633 forAll(newFlipMode, i)
1635 newFlipMode[i] = flipMode[zoneI][newToOld[i]];
1637 flipMode[zoneI].transfer(newFlipMode);
1641 // So now we both have old zones and the new addressing.
1642 // Invert the addressing to get faceZoneFaceMap.
1643 forAll(addressing, zoneI)
1645 const faceZone& oldZone = faceZones[zoneI];
1646 const labelList& newZoneAddr = addressing[zoneI];
1648 labelList& curFzFaceRnb = faceZoneFaceMap[zoneI];
1650 curFzFaceRnb.setSize(newZoneAddr.size());
1652 forAll(newZoneAddr, i)
1654 if (newZoneAddr[i] < faceMap_.size())
1657 oldZone.whichFace(faceMap_[newZoneAddr[i]]);
1661 curFzFaceRnb[i] = -1;
1667 // Reset the addresing on the zone
1668 newMesh.faceZones().clearAddressing();
1669 forAll(newMesh.faceZones(), zoneI)
1673 Pout<< "faceZone:" << zoneI
1674 << " name:" << newMesh.faceZones()[zoneI].name()
1675 << " size:" << addressing[zoneI].size()
1679 newMesh.faceZones()[zoneI].resetAddressing
1691 cellZoneMap.setSize(mesh.cellZones().size());
1693 const cellZoneMesh& cellZones = mesh.cellZones();
1695 labelList nCells(cellZones.size(), 0);
1697 forAll(cellZone_, cellI)
1699 label zoneI = cellZone_[cellI];
1701 if (zoneI >= cellZones.size())
1705 "resetZones(const polyMesh&, polyMesh&, labelListList&"
1706 "labelListList&, labelListList&)"
1707 ) << "Illegal zoneID " << zoneI << " for cell "
1708 << cellI << abort(FatalError);
1717 labelListList addressing(cellZones.size());
1718 forAll(addressing, zoneI)
1720 addressing[zoneI].setSize(nCells[zoneI]);
1724 forAll(cellZone_, cellI)
1726 label zoneI = cellZone_[cellI];
1730 addressing[zoneI][nCells[zoneI]++] = cellI;
1733 // Sort the addressing
1734 forAll(addressing, zoneI)
1736 stableSort(addressing[zoneI]);
1739 // So now we both have old zones and the new addressing.
1740 // Invert the addressing to get cellZoneMap.
1741 forAll(addressing, zoneI)
1743 const cellZone& oldZone = cellZones[zoneI];
1744 const labelList& newZoneAddr = addressing[zoneI];
1746 labelList& curCellRnb = cellZoneMap[zoneI];
1748 curCellRnb.setSize(newZoneAddr.size());
1750 forAll(newZoneAddr, i)
1752 if (newZoneAddr[i] < cellMap_.size())
1755 oldZone.whichCell(cellMap_[newZoneAddr[i]]);
1764 // Reset the addresing on the zone
1765 newMesh.cellZones().clearAddressing();
1766 forAll(newMesh.cellZones(), zoneI)
1770 Pout<< "cellZone:" << zoneI
1771 << " name:" << newMesh.cellZones()[zoneI].name()
1772 << " size:" << addressing[zoneI].size()
1776 newMesh.cellZones()[zoneI] = addressing[zoneI];
1782 void Foam::polyTopoChange::calcFaceZonePointMap
1784 const polyMesh& mesh,
1785 const List<Map<label> >& oldFaceZoneMeshPointMaps,
1786 labelListList& faceZonePointMap
1789 const faceZoneMesh& faceZones = mesh.faceZones();
1791 faceZonePointMap.setSize(faceZones.size());
1793 forAll(faceZones, zoneI)
1795 const faceZone& newZone = faceZones[zoneI];
1797 const labelList& newZoneMeshPoints = newZone().meshPoints();
1799 const Map<label>& oldZoneMeshPointMap = oldFaceZoneMeshPointMaps[zoneI];
1801 labelList& curFzPointRnb = faceZonePointMap[zoneI];
1803 curFzPointRnb.setSize(newZoneMeshPoints.size());
1805 forAll(newZoneMeshPoints, pointI)
1807 if (newZoneMeshPoints[pointI] < pointMap_.size())
1809 Map<label>::const_iterator ozmpmIter =
1810 oldZoneMeshPointMap.find
1812 pointMap_[newZoneMeshPoints[pointI]]
1815 if (ozmpmIter != oldZoneMeshPointMap.end())
1817 curFzPointRnb[pointI] = ozmpmIter();
1821 curFzPointRnb[pointI] = -1;
1826 curFzPointRnb[pointI] = -1;
1833 Foam::face Foam::polyTopoChange::rotateFace
1839 face newF(f.size());
1843 label fp1 = (fp + nPos) % f.size();
1857 void Foam::polyTopoChange::reorderCoupledFaces
1859 const bool syncParallel,
1860 const polyBoundaryMesh& boundary,
1861 const labelList& patchStarts,
1862 const labelList& patchSizes,
1863 const pointField& points
1866 // Mapping for faces (old to new). Extends over all mesh faces for
1867 // convenience (could be just the external faces)
1868 labelList oldToNew(identity(faces_.size()));
1870 // Rotation on new faces.
1871 labelList rotation(faces_.size(), 0);
1874 forAll(boundary, patchI)
1876 if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1878 boundary[patchI].initOrder
1894 // Receive and calculate ordering
1896 bool anyChanged = false;
1898 forAll(boundary, patchI)
1900 if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1902 labelList patchFaceMap(patchSizes[patchI], -1);
1903 labelList patchFaceRotation(patchSizes[patchI], 0);
1905 bool changed = boundary[patchI].order
1923 // Merge patch face reordering into mesh face reordering table
1924 label start = patchStarts[patchI];
1926 forAll(patchFaceMap, patchFaceI)
1928 oldToNew[patchFaceI + start] =
1929 start + patchFaceMap[patchFaceI];
1932 forAll(patchFaceRotation, patchFaceI)
1934 rotation[patchFaceI + start] =
1935 patchFaceRotation[patchFaceI];
1945 reduce(anyChanged, orOp<bool>());
1950 // Reorder faces according to oldToNew.
1951 reorderCompactFaces(oldToNew.size(), oldToNew);
1953 // Rotate faces (rotation is already in new face indices).
1954 forAll(rotation, faceI)
1956 if (rotation[faceI] != 0)
1958 faces_[faceI] = rotateFace(faces_[faceI], rotation[faceI]);
1965 void Foam::polyTopoChange::compactAndReorder
1967 const polyMesh& mesh,
1968 const bool syncParallel,
1969 const bool orderCells,
1970 const bool orderPoints,
1972 label& nInternalPoints,
1973 pointField& newPoints,
1974 labelList& patchSizes,
1975 labelList& patchStarts,
1976 List<objectMap>& pointsFromPoints,
1977 List<objectMap>& facesFromPoints,
1978 List<objectMap>& facesFromEdges,
1979 List<objectMap>& facesFromFaces,
1980 List<objectMap>& cellsFromPoints,
1981 List<objectMap>& cellsFromEdges,
1982 List<objectMap>& cellsFromFaces,
1983 List<objectMap>& cellsFromCells,
1984 List<Map<label> >& oldPatchMeshPointMaps,
1985 labelList& oldPatchNMeshPoints,
1986 labelList& oldPatchStarts,
1987 List<Map<label> >& oldFaceZoneMeshPointMaps
1990 if (mesh.boundaryMesh().size() != nPatches_)
1992 FatalErrorIn("polyTopoChange::compactAndReorder(..)")
1993 << "polyTopoChange was constructed with a mesh with "
1994 << nPatches_ << " patches." << endl
1995 << "The mesh now provided has a different number of patches "
1996 << mesh.boundaryMesh().size()
1997 << " which is illegal" << endl
1998 << abort(FatalError);
2001 // Remove any holes from points/faces/cells and sort faces.
2002 // Sets nActiveFaces_.
2003 compact(orderCells, orderPoints, nInternalPoints, patchSizes, patchStarts);
2005 // Transfer points to pointField. points_ are now cleared!
2006 // Only done since e.g. reorderCoupledFaces requires pointField.
2007 newPoints.transfer(points_);
2009 // Reorder any coupled faces
2013 mesh.boundaryMesh(),
2020 // Calculate inflation/merging maps
2021 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2022 // These are for the new face(/point/cell) the old faces whose value
2024 // averaged/summed to get the new value. These old faces come either from
2025 // merged old faces (face remove into other face),
2026 // the old edgeFaces (inflate from edge) or the old pointFaces (inflate
2027 // from point). As an additional complexity will use only internal faces
2028 // to create new value for internal face and vice versa only patch
2029 // faces to to create patch face value.
2031 // For point only point merging
2039 calcFaceInflationMaps
2047 calcCellInflationMaps
2056 // Clear inflation info
2058 faceFromPoint_.clearStorage();
2059 faceFromEdge_.clearStorage();
2061 cellFromPoint_.clearStorage();
2062 cellFromEdge_.clearStorage();
2063 cellFromFace_.clearStorage();
2067 const polyBoundaryMesh& boundary = mesh.boundaryMesh();
2069 // Grab patch mesh point maps
2070 oldPatchMeshPointMaps.setSize(boundary.size());
2071 oldPatchNMeshPoints.setSize(boundary.size());
2072 oldPatchStarts.setSize(boundary.size());
2074 forAll(boundary, patchI)
2076 // Copy old face zone mesh point maps
2077 oldPatchMeshPointMaps[patchI] = boundary[patchI].meshPointMap();
2078 oldPatchNMeshPoints[patchI] = boundary[patchI].meshPoints().size();
2079 oldPatchStarts[patchI] = boundary[patchI].start();
2082 // Grab old face zone mesh point maps.
2083 // These need to be saved before resetting the mesh and are used
2084 // later on to calculate the faceZone pointMaps.
2085 oldFaceZoneMeshPointMaps.setSize(mesh.faceZones().size());
2087 forAll(mesh.faceZones(), zoneI)
2089 const faceZone& oldZone = mesh.faceZones()[zoneI];
2091 oldFaceZoneMeshPointMaps[zoneI] = oldZone().meshPointMap();
2096 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
2098 // Construct from components
2099 Foam::polyTopoChange::polyTopoChange(const label nPatches, const bool strict)
2102 nPatches_(nPatches),
2105 reversePointMap_(0),
2129 // Construct from components
2130 Foam::polyTopoChange::polyTopoChange
2132 const polyMesh& mesh,
2140 reversePointMap_(0),
2165 identity(mesh.boundaryMesh().size()),
2166 identity(mesh.pointZones().size()),
2167 identity(mesh.faceZones().size()),
2168 identity(mesh.cellZones().size())
2173 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2175 void Foam::polyTopoChange::clear()
2177 points_.clearStorage();
2178 pointMap_.clearStorage();
2179 reversePointMap_.clearStorage();
2180 pointZone_.clearStorage();
2181 retiredPoints_.clearStorage();
2183 faces_.clearStorage();
2184 region_.clearStorage();
2185 faceOwner_.clearStorage();
2186 faceNeighbour_.clearStorage();
2187 faceMap_.clearStorage();
2188 reverseFaceMap_.clearStorage();
2189 faceFromPoint_.clearStorage();
2190 faceFromEdge_.clearStorage();
2191 flipFaceFlux_.clearStorage();
2192 faceZone_.clearStorage();
2193 faceZoneFlip_.clearStorage();
2196 cellMap_.clearStorage();
2197 reverseCellMap_.clearStorage();
2198 cellZone_.clearStorage();
2199 cellFromPoint_.clearStorage();
2200 cellFromEdge_.clearStorage();
2201 cellFromFace_.clearStorage();
2205 void Foam::polyTopoChange::addMesh
2207 const polyMesh& mesh,
2208 const labelList& patchMap,
2209 const labelList& pointZoneMap,
2210 const labelList& faceZoneMap,
2211 const labelList& cellZoneMap
2214 label maxRegion = nPatches_ - 1;
2217 maxRegion = max(maxRegion, patchMap[i]);
2219 nPatches_ = maxRegion + 1;
2224 const pointField& points = mesh.points();
2225 const pointZoneMesh& pointZones = mesh.pointZones();
2228 points_.setCapacity(points_.size() + points.size());
2229 pointMap_.setCapacity(pointMap_.size() + points.size());
2230 reversePointMap_.setCapacity(reversePointMap_.size() + points.size());
2231 pointZone_.setCapacity(pointZone_.size() + points.size());
2233 // Precalc offset zones
2234 labelList newZoneID(points.size(), -1);
2236 forAll(pointZones, zoneI)
2238 const labelList& pointLabels = pointZones[zoneI];
2240 forAll(pointLabels, j)
2242 newZoneID[pointLabels[j]] = pointZoneMap[zoneI];
2246 // Add points in mesh order
2247 for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
2261 const cellZoneMesh& cellZones = mesh.cellZones();
2265 // Note: polyMesh does not allow retired cells anymore. So allCells
2266 // always equals nCells
2267 label nAllCells = mesh.nCells();
2269 cellMap_.setCapacity(cellMap_.size() + nAllCells);
2270 reverseCellMap_.setCapacity(reverseCellMap_.size() + nAllCells);
2271 cellFromPoint_.resize(cellFromPoint_.size() + nAllCells/100);
2272 cellFromEdge_.resize(cellFromEdge_.size() + nAllCells/100);
2273 cellFromFace_.resize(cellFromFace_.size() + nAllCells/100);
2274 cellZone_.setCapacity(cellZone_.size() + nAllCells);
2277 // Precalc offset zones
2278 labelList newZoneID(nAllCells, -1);
2280 forAll(cellZones, zoneI)
2282 const labelList& cellLabels = cellZones[zoneI];
2284 forAll(cellLabels, j)
2286 label cellI = cellLabels[j];
2288 if (newZoneID[cellI] != -1)
2292 "polyTopoChange::addMesh"
2293 "(const polyMesh&, const labelList&,"
2294 "const labelList&, const labelList&,"
2296 ) << "Cell:" << cellI
2297 << " centre:" << mesh.cellCentres()[cellI]
2298 << " is in two zones:"
2299 << cellZones[newZoneID[cellI]].name()
2300 << " and " << cellZones[zoneI].name() << endl
2301 << " This is not supported."
2302 << " Continuing with first zone only." << endl;
2306 newZoneID[cellI] = cellZoneMap[zoneI];
2311 // Add cells in mesh order
2312 for (label cellI = 0; cellI < nAllCells; cellI++)
2314 // Add cell from cell
2315 addCell(-1, -1, -1, cellI, newZoneID[cellI]);
2321 const polyBoundaryMesh& patches = mesh.boundaryMesh();
2322 const faceList& faces = mesh.faces();
2323 const labelList& faceOwner = mesh.faceOwner();
2324 const labelList& faceNeighbour = mesh.faceNeighbour();
2325 const faceZoneMesh& faceZones = mesh.faceZones();
2328 label nAllFaces = mesh.faces().size();
2330 faces_.setCapacity(faces_.size() + nAllFaces);
2331 region_.setCapacity(region_.size() + nAllFaces);
2332 faceOwner_.setCapacity(faceOwner_.size() + nAllFaces);
2333 faceNeighbour_.setCapacity(faceNeighbour_.size() + nAllFaces);
2334 faceMap_.setCapacity(faceMap_.size() + nAllFaces);
2335 reverseFaceMap_.setCapacity(reverseFaceMap_.size() + nAllFaces);
2336 faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/100);
2337 faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/100);
2338 flipFaceFlux_.setCapacity(flipFaceFlux_.size() + nAllFaces);
2339 faceZone_.setCapacity(faceZone_.size() + nAllFaces);
2340 faceZoneFlip_.setCapacity(faceZoneFlip_.size() + nAllFaces);
2343 // Precalc offset zones
2344 labelList newZoneID(nAllFaces, -1);
2345 boolList zoneFlip(nAllFaces, false);
2347 forAll(faceZones, zoneI)
2349 const labelList& faceLabels = faceZones[zoneI];
2350 const boolList& flipMap = faceZones[zoneI].flipMap();
2352 forAll(faceLabels, j)
2354 newZoneID[faceLabels[j]] = faceZoneMap[zoneI];
2355 zoneFlip[faceLabels[j]] = flipMap[j];
2359 // Add faces in mesh order
2361 // 1. Internal faces
2362 for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
2368 faceNeighbour[faceI],
2369 -1, // masterPointID
2371 faceI, // masterFaceID
2372 false, // flipFaceFlux
2374 newZoneID[faceI], // zoneID
2375 zoneFlip[faceI] // zoneFlip
2380 forAll(patches, patchI)
2382 const polyPatch& pp = patches[patchI];
2384 if (pp.start() != faces_.size())
2388 "polyTopoChange::polyTopoChange"
2389 "(const polyMesh& mesh, const bool strict)"
2391 << "Patch " << pp.name() << " starts at " << pp.start()
2393 << "Current face counter at " << faces_.size() << endl
2394 << "Are patches in incremental order?"
2395 << abort(FatalError);
2397 forAll(pp, patchFaceI)
2399 label faceI = pp.start() + patchFaceI;
2406 -1, // masterPointID
2408 faceI, // masterFaceID
2409 false, // flipFaceFlux
2410 patchMap[patchI], // patchID
2411 newZoneID[faceI], // zoneID
2412 zoneFlip[faceI] // zoneFlip
2420 void Foam::polyTopoChange::setCapacity
2422 const label nPoints,
2427 points_.setCapacity(nPoints);
2428 pointMap_.setCapacity(nPoints);
2429 reversePointMap_.setCapacity(nPoints);
2430 pointZone_.setCapacity(nPoints);
2432 faces_.setCapacity(nFaces);
2433 region_.setCapacity(nFaces);
2434 faceOwner_.setCapacity(nFaces);
2435 faceNeighbour_.setCapacity(nFaces);
2436 faceMap_.setCapacity(nFaces);
2437 reverseFaceMap_.setCapacity(nFaces);
2438 faceFromPoint_.resize(faceFromPoint_.size() + nFaces/100);
2439 faceFromEdge_.resize(faceFromEdge_.size() + nFaces/100);
2440 flipFaceFlux_.setCapacity(nFaces);
2441 faceZone_.setCapacity(nFaces);
2442 faceZoneFlip_.setCapacity(nFaces);
2444 cellMap_.setCapacity(nCells);
2445 reverseCellMap_.setCapacity(nCells);
2446 cellFromPoint_.resize(cellFromPoint_.size() + nCells/100);
2447 cellFromEdge_.resize(cellFromEdge_.size() + nCells/100);
2448 cellFromFace_.resize(cellFromFace_.size() + nCells/100);
2449 cellZone_.setCapacity(nCells);
2453 Foam::label Foam::polyTopoChange::setAction(const topoAction& action)
2455 if (isType<polyAddPoint>(action))
2457 const polyAddPoint& pap = refCast<const polyAddPoint>(action);
2462 pap.masterPointID(),
2467 else if (isType<polyModifyPoint>(action))
2469 const polyModifyPoint& pmp = refCast<const polyModifyPoint>(action);
2481 else if (isType<polyRemovePoint>(action))
2483 const polyRemovePoint& prp = refCast<const polyRemovePoint>(action);
2485 removePoint(prp.pointID(), prp.mergePointID());
2489 else if (isType<polyAddFace>(action))
2491 const polyAddFace& paf = refCast<const polyAddFace>(action);
2498 paf.masterPointID(),
2507 else if (isType<polyModifyFace>(action))
2509 const polyModifyFace& pmf = refCast<const polyModifyFace>(action);
2525 else if (isType<polyRemoveFace>(action))
2527 const polyRemoveFace& prf = refCast<const polyRemoveFace>(action);
2529 removeFace(prf.faceID(), prf.mergeFaceID());
2533 else if (isType<polyAddCell>(action))
2535 const polyAddCell& pac = refCast<const polyAddCell>(action);
2539 pac.masterPointID(),
2546 else if (isType<polyModifyCell>(action))
2548 const polyModifyCell& pmc = refCast<const polyModifyCell>(action);
2550 if (pmc.removeFromZone())
2552 modifyCell(pmc.cellID(), -1);
2556 modifyCell(pmc.cellID(), pmc.zoneID());
2561 else if (isType<polyRemoveCell>(action))
2563 const polyRemoveCell& prc = refCast<const polyRemoveCell>(action);
2565 removeCell(prc.cellID(), prc.mergeCellID());
2573 "label polyTopoChange::setAction(const topoAction& action)"
2574 ) << "Unknown type of topoChange: " << action.type()
2575 << abort(FatalError);
2577 // Dummy return to keep compiler happy
2583 Foam::label Foam::polyTopoChange::addPoint
2586 const label masterPointID,
2591 label pointI = points_.size();
2594 pointMap_.append(masterPointID);
2595 reversePointMap_.append(pointI);
2596 pointZone_.append(zoneID);
2600 retiredPoints_.insert(pointI);
2607 void Foam::polyTopoChange::modifyPoint
2611 const label newZoneID,
2615 if (pointI < 0 || pointI >= points_.size())
2619 "polyTopoChange::modifyPoint(const label, const point&)"
2620 ) << "illegal point label " << pointI << endl
2621 << "Valid point labels are 0 .. " << points_.size()-1
2622 << abort(FatalError);
2624 if (pointRemoved(pointI) || pointMap_[pointI] == -1)
2628 "polyTopoChange::modifyPoint(const label, const point&)"
2629 ) << "point " << pointI << " already marked for removal"
2630 << abort(FatalError);
2632 points_[pointI] = pt;
2633 pointZone_[pointI] = newZoneID;
2637 retiredPoints_.erase(pointI);
2641 retiredPoints_.insert(pointI);
2646 void Foam::polyTopoChange::movePoints(const pointField& newPoints)
2648 if (newPoints.size() != points_.size())
2650 FatalErrorIn("polyTopoChange::movePoints(const pointField&)")
2651 << "illegal pointField size." << endl
2652 << "Size:" << newPoints.size() << endl
2653 << "Points in mesh:" << points_.size()
2654 << abort(FatalError);
2657 forAll(points_, pointI)
2659 points_[pointI] = newPoints[pointI];
2664 void Foam::polyTopoChange::removePoint
2667 const label mergePointI
2670 if (pointI < 0 || pointI >= points_.size())
2672 FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2673 << "illegal point label " << pointI << endl
2674 << "Valid point labels are 0 .. " << points_.size()-1
2675 << abort(FatalError);
2681 && (pointRemoved(pointI) || pointMap_[pointI] == -1)
2684 FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2685 << "point " << pointI << " already marked for removal" << nl
2686 << "Point:" << points_[pointI] << " pointMap:" << pointMap_[pointI]
2687 << abort(FatalError);
2690 if (pointI == mergePointI)
2692 FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2693 << "Cannot remove/merge point " << pointI << " onto itself."
2694 << abort(FatalError);
2697 points_[pointI] = greatPoint;
2698 pointMap_[pointI] = -1;
2699 if (mergePointI >= 0)
2701 reversePointMap_[pointI] = -mergePointI-2;
2705 reversePointMap_[pointI] = -1;
2707 pointZone_[pointI] = -1;
2708 retiredPoints_.erase(pointI);
2712 Foam::label Foam::polyTopoChange::addFace
2717 const label masterPointID,
2718 const label masterEdgeID,
2719 const label masterFaceID,
2720 const bool flipFaceFlux,
2721 const label patchID,
2729 checkFace(f, -1, own, nei, patchID, zoneID);
2732 label faceI = faces_.size();
2735 region_.append(patchID);
2736 faceOwner_.append(own);
2737 faceNeighbour_.append(nei);
2739 if (masterPointID >= 0)
2741 faceMap_.append(-1);
2742 faceFromPoint_.insert(faceI, masterPointID);
2744 else if (masterEdgeID >= 0)
2746 faceMap_.append(-1);
2747 faceFromEdge_.insert(faceI, masterEdgeID);
2749 else if (masterFaceID >= 0)
2751 faceMap_.append(masterFaceID);
2755 // Allow inflate-from-nothing?
2756 //FatalErrorIn("polyTopoChange::addFace")
2757 // << "Need to specify a master point, edge or face"
2758 // << "face:" << f << " own:" << own << " nei:" << nei
2759 // << abort(FatalError);
2760 faceMap_.append(-1);
2762 reverseFaceMap_.append(faceI);
2764 flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2765 faceZone_.append(zoneID);
2766 faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2772 void Foam::polyTopoChange::modifyFace
2778 const bool flipFaceFlux,
2779 const label patchID,
2787 checkFace(f, faceI, own, nei, patchID, zoneID);
2791 faceOwner_[faceI] = own;
2792 faceNeighbour_[faceI] = nei;
2793 region_[faceI] = patchID;
2795 flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2796 faceZone_[faceI] = zoneID;
2797 faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2801 void Foam::polyTopoChange::removeFace(const label faceI, const label mergeFaceI)
2803 if (faceI < 0 || faceI >= faces_.size())
2805 FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2806 << "illegal face label " << faceI << endl
2807 << "Valid face labels are 0 .. " << faces_.size()-1
2808 << abort(FatalError);
2814 && (faceRemoved(faceI) || faceMap_[faceI] == -1)
2817 FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2819 << " already marked for removal"
2820 << abort(FatalError);
2823 faces_[faceI].setSize(0);
2824 region_[faceI] = -1;
2825 faceOwner_[faceI] = -1;
2826 faceNeighbour_[faceI] = -1;
2827 faceMap_[faceI] = -1;
2828 if (mergeFaceI >= 0)
2830 reverseFaceMap_[faceI] = -mergeFaceI-2;
2834 reverseFaceMap_[faceI] = -1;
2836 faceFromEdge_.erase(faceI);
2837 faceFromPoint_.erase(faceI);
2838 flipFaceFlux_[faceI] = 0;
2839 faceZone_[faceI] = -1;
2840 faceZoneFlip_[faceI] = 0;
2844 Foam::label Foam::polyTopoChange::addCell
2846 const label masterPointID,
2847 const label masterEdgeID,
2848 const label masterFaceID,
2849 const label masterCellID,
2853 label cellI = cellMap_.size();
2855 if (masterPointID >= 0)
2857 cellMap_.append(-1);
2858 cellFromPoint_.insert(cellI, masterPointID);
2860 else if (masterEdgeID >= 0)
2862 cellMap_.append(-1);
2863 cellFromEdge_.insert(cellI, masterEdgeID);
2865 else if (masterFaceID >= 0)
2867 cellMap_.append(-1);
2868 cellFromFace_.insert(cellI, masterFaceID);
2872 cellMap_.append(masterCellID);
2874 reverseCellMap_.append(cellI);
2875 cellZone_.append(zoneID);
2881 void Foam::polyTopoChange::modifyCell
2887 cellZone_[cellI] = zoneID;
2891 void Foam::polyTopoChange::removeCell(const label cellI, const label mergeCellI)
2893 if (cellI < 0 || cellI >= cellMap_.size())
2895 FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2896 << "illegal cell label " << cellI << endl
2897 << "Valid cell labels are 0 .. " << cellMap_.size()-1
2898 << abort(FatalError);
2901 if (strict_ && cellMap_[cellI] == -2)
2903 FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2905 << " already marked for removal"
2906 << abort(FatalError);
2909 cellMap_[cellI] = -2;
2910 if (mergeCellI >= 0)
2912 reverseCellMap_[cellI] = -mergeCellI-2;
2916 reverseCellMap_[cellI] = -1;
2918 cellFromPoint_.erase(cellI);
2919 cellFromEdge_.erase(cellI);
2920 cellFromFace_.erase(cellI);
2921 cellZone_[cellI] = -1;
2925 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
2929 const bool syncParallel,
2930 const bool orderCells,
2931 const bool orderPoints
2936 Pout<< "polyTopoChange::changeMesh"
2937 << "(polyMesh&, const bool, const bool, const bool, const bool)"
2943 Pout<< "Old mesh:" << nl;
2944 writeMeshStats(mesh, Pout);
2948 pointField newPoints;
2949 // number of internal points
2950 label nInternalPoints;
2952 labelList patchSizes;
2953 labelList patchStarts;
2955 List<objectMap> pointsFromPoints;
2956 List<objectMap> facesFromPoints;
2957 List<objectMap> facesFromEdges;
2958 List<objectMap> facesFromFaces;
2959 List<objectMap> cellsFromPoints;
2960 List<objectMap> cellsFromEdges;
2961 List<objectMap> cellsFromFaces;
2962 List<objectMap> cellsFromCells;
2964 List<Map<label> > oldPatchMeshPointMaps;
2965 labelList oldPatchNMeshPoints;
2966 labelList oldPatchStarts;
2967 List<Map<label> > oldFaceZoneMeshPointMaps;
2969 // Compact, reorder patch faces and calculate mesh/patch maps.
2989 oldPatchMeshPointMaps,
2990 oldPatchNMeshPoints,
2992 oldFaceZoneMeshPointMaps
2995 const label nOldPoints(mesh.nPoints());
2996 const label nOldFaces(mesh.nFaces());
2997 const label nOldCells(mesh.nCells());
3002 // This will invalidate any addressing so better make sure you have
3003 // all the information you need!!!
3007 // Keep (renumbered) mesh points, store new points in map for inflation
3008 // (appended points (i.e. from nowhere) get value zero)
3009 pointField renumberedMeshPoints(newPoints.size());
3011 forAll(pointMap_, newPointI)
3013 label oldPointI = pointMap_[newPointI];
3017 renumberedMeshPoints[newPointI] = mesh.points()[oldPointI];
3021 renumberedMeshPoints[newPointI] = vector::zero;
3025 mesh.resetPrimitives
3027 xferMove(renumberedMeshPoints),
3030 faceNeighbour_.xfer(),
3036 mesh.changing(true);
3041 mesh.resetPrimitives
3043 xferMove(newPoints),
3046 faceNeighbour_.xfer(),
3051 mesh.changing(true);
3054 // Clear out primitives
3056 retiredPoints_.clearStorage();
3057 region_.clearStorage();
3063 // Some stats on changes
3064 label nAdd, nInflate, nMerge, nRemove;
3065 countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3067 << " added(from point):" << nAdd
3068 << " added(from nothing):" << nInflate
3069 << " merged(into other point):" << nMerge
3070 << " removed:" << nRemove
3073 countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3075 << " added(from face):" << nAdd
3076 << " added(inflated):" << nInflate
3077 << " merged(into other face):" << nMerge
3078 << " removed:" << nRemove
3081 countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3083 << " added(from cell):" << nAdd
3084 << " added(inflated):" << nInflate
3085 << " merged(into other cell):" << nMerge
3086 << " removed:" << nRemove
3093 Pout<< "New mesh:" << nl;
3094 writeMeshStats(mesh, Pout);
3101 // Inverse of point/face/cell zone addressing.
3102 // For every preserved point/face/cells in zone give the old position.
3103 // For added points, the index is set to -1
3104 labelListList pointZoneMap(mesh.pointZones().size());
3105 labelListList faceZoneFaceMap(mesh.faceZones().size());
3106 labelListList cellZoneMap(mesh.cellZones().size());
3108 resetZones(mesh, mesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3112 pointZone_.clearStorage();
3113 faceZone_.clearStorage();
3114 faceZoneFlip_.clearStorage();
3115 cellZone_.clearStorage();
3119 // Patch point renumbering
3120 // For every preserved point on a patch give the old position.
3121 // For added points, the index is set to -1
3122 labelListList patchPointMap(mesh.boundaryMesh().size());
3125 oldPatchMeshPointMaps,
3126 mesh.boundaryMesh(),
3130 // Create the face zone mesh point renumbering
3131 labelListList faceZonePointMap(mesh.faceZones().size());
3132 calcFaceZonePointMap(mesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3134 labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3136 return autoPtr<mapPolyMesh>
3173 newPoints, // if empty signals no inflation.
3175 oldPatchNMeshPoints,
3176 true // steal storage.
3180 // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3185 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
3187 autoPtr<fvMesh>& newMeshPtr,
3190 const bool syncParallel,
3191 const bool orderCells,
3192 const bool orderPoints
3197 Pout<< "polyTopoChange::changeMesh"
3198 << "(autoPtr<fvMesh>&, const IOobject&, const fvMesh&"
3199 << ", const bool, const bool, const bool)"
3205 Pout<< "Old mesh:" << nl;
3206 writeMeshStats(mesh, Pout);
3210 pointField newPoints;
3211 // number of internal points
3212 label nInternalPoints;
3214 labelList patchSizes;
3215 labelList patchStarts;
3217 List<objectMap> pointsFromPoints;
3218 List<objectMap> facesFromPoints;
3219 List<objectMap> facesFromEdges;
3220 List<objectMap> facesFromFaces;
3221 List<objectMap> cellsFromPoints;
3222 List<objectMap> cellsFromEdges;
3223 List<objectMap> cellsFromFaces;
3224 List<objectMap> cellsFromCells;
3227 List<Map<label> > oldPatchMeshPointMaps;
3228 labelList oldPatchNMeshPoints;
3229 labelList oldPatchStarts;
3230 List<Map<label> > oldFaceZoneMeshPointMaps;
3232 // Compact, reorder patch faces and calculate mesh/patch maps.
3252 oldPatchMeshPointMaps,
3253 oldPatchNMeshPoints,
3255 oldFaceZoneMeshPointMaps
3258 const label nOldPoints(mesh.nPoints());
3259 const label nOldFaces(mesh.nFaces());
3260 const label nOldCells(mesh.nCells());
3271 xferMove(newPoints),
3274 faceNeighbour_.xfer()
3277 fvMesh& newMesh = newMeshPtr();
3279 // Clear out primitives
3281 retiredPoints_.clearStorage();
3282 region_.clearStorage();
3288 // Some stats on changes
3289 label nAdd, nInflate, nMerge, nRemove;
3290 countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3292 << " added(from point):" << nAdd
3293 << " added(from nothing):" << nInflate
3294 << " merged(into other point):" << nMerge
3295 << " removed:" << nRemove
3298 countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3300 << " added(from face):" << nAdd
3301 << " added(inflated):" << nInflate
3302 << " merged(into other face):" << nMerge
3303 << " removed:" << nRemove
3306 countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3308 << " added(from cell):" << nAdd
3309 << " added(inflated):" << nInflate
3310 << " merged(into other cell):" << nMerge
3311 << " removed:" << nRemove
3318 const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
3320 List<polyPatch*> newBoundary(oldPatches.size());
3322 forAll(oldPatches, patchI)
3324 newBoundary[patchI] = oldPatches[patchI].clone
3326 newMesh.boundaryMesh(),
3332 newMesh.addFvPatches(newBoundary);
3339 // Start off from empty zones.
3340 const pointZoneMesh& oldPointZones = mesh.pointZones();
3341 List<pointZone*> pZonePtrs(oldPointZones.size());
3343 forAll(oldPointZones, i)
3345 pZonePtrs[i] = new pointZone
3347 oldPointZones[i].name(),
3350 newMesh.pointZones()
3355 const faceZoneMesh& oldFaceZones = mesh.faceZones();
3356 List<faceZone*> fZonePtrs(oldFaceZones.size());
3358 forAll(oldFaceZones, i)
3360 fZonePtrs[i] = new faceZone
3362 oldFaceZones[i].name(),
3371 const cellZoneMesh& oldCellZones = mesh.cellZones();
3372 List<cellZone*> cZonePtrs(oldCellZones.size());
3374 forAll(oldCellZones, i)
3376 cZonePtrs[i] = new cellZone
3378 oldCellZones[i].name(),
3386 newMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
3388 // Inverse of point/face/cell zone addressing.
3389 // For every preserved point/face/cells in zone give the old position.
3390 // For added points, the index is set to -1
3391 labelListList pointZoneMap(mesh.pointZones().size());
3392 labelListList faceZoneFaceMap(mesh.faceZones().size());
3393 labelListList cellZoneMap(mesh.cellZones().size());
3395 resetZones(mesh, newMesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3399 pointZone_.clearStorage();
3400 faceZone_.clearStorage();
3401 faceZoneFlip_.clearStorage();
3402 cellZone_.clearStorage();
3405 // Patch point renumbering
3406 // For every preserved point on a patch give the old position.
3407 // For added points, the index is set to -1
3408 labelListList patchPointMap(newMesh.boundaryMesh().size());
3411 oldPatchMeshPointMaps,
3412 newMesh.boundaryMesh(),
3416 // Create the face zone mesh point renumbering
3417 labelListList faceZonePointMap(newMesh.faceZones().size());
3418 calcFaceZonePointMap(newMesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3422 Pout<< "New mesh:" << nl;
3423 writeMeshStats(mesh, Pout);
3426 labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3428 return autoPtr<mapPolyMesh>
3465 newPoints, // if empty signals no inflation.
3467 oldPatchNMeshPoints,
3468 true // steal storage.
3472 // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3477 // ************************************************************************* //