61 void Foam::polyTopoChange::renumberReverseMap
64 DynamicList<label>& elems
69 label val = elems[elemI];
73 elems[elemI] = map[val];
77 label mergedVal = -val-2;
78 elems[elemI] = -map[mergedVal]-2;
94 label newElem = map[iter.key()];
98 newElems.insert(newElem);
102 elems.transfer(newElems);
107 void Foam::polyTopoChange::renumberCompact
117 label newVal = map[elems[elemI]];
121 elems[newElemI++] = newVal;
124 elems.setSize(newElemI);
128 void Foam::polyTopoChange::countMap
145 label oldCellI = map[newCellI];
149 if (reverseMap[oldCellI] == newCellI)
159 else if (oldCellI == -1)
171 forAll(reverseMap, oldCellI)
173 label newCellI = reverseMap[oldCellI];
179 else if (newCellI == -1)
210 void Foam::polyTopoChange::writeMeshStats(
const polyMesh&
mesh, Ostream& os)
212 const polyBoundaryMesh&
patches = mesh.boundaryMesh();
218 patchSizes[patchI] = patches[patchI].size();
219 patchStarts[patchI] = patches[patchI].start();
222 os <<
" Points : " << mesh.nPoints() <<
nl
223 <<
" Faces : " << mesh.nFaces() <<
nl
224 <<
" Cells : " << mesh.nCells() <<
nl
225 <<
" PatchSizes : " << patchSizes <<
nl
226 <<
" PatchStarts : " << patchStarts <<
nl
231 void Foam::polyTopoChange::getMergeSets
235 List<objectMap>& cellsFromCells
241 forAll(reverseCellMap, oldCellI)
243 label newCellI = reverseCellMap[oldCellI];
247 label mergeCellI = -newCellI-2;
249 nMerged[mergeCellI]++;
254 labelList cellToMergeSet(cellMap.size(), -1);
260 if (nMerged[cellI] > 1)
262 cellToMergeSet[cellI] = nSets++;
272 cellsFromCells.setSize(nSets);
274 forAll(reverseCellMap, oldCellI)
276 label newCellI = reverseCellMap[oldCellI];
280 label mergeCellI = -newCellI-2;
284 label setI = cellToMergeSet[mergeCellI];
286 objectMap& mergeSet = cellsFromCells[setI];
288 if (mergeSet.masterObjects().empty())
292 mergeSet.index() = mergeCellI;
293 mergeSet.masterObjects().setSize(nMerged[mergeCellI]);
296 mergeSet.masterObjects()[0] = cellMap[mergeCellI];
299 mergeSet.masterObjects()[1] = oldCellI;
301 nMerged[mergeCellI] = 2;
305 mergeSet.masterObjects()[nMerged[mergeCellI]++] = oldCellI;
312 void Foam::polyTopoChange::checkFace
324 if (own == -1 && zoneI != -1)
328 else if (patchI == -1 || patchI >= nPatches_)
332 "polyTopoChange::checkFace(const face&, const label"
333 ", const label, const label, const label)"
334 ) <<
"Face has no neighbour (so external) but does not have"
335 <<
" a valid patch" <<
nl
337 <<
" faceI(-1 if added face):" << faceI
338 <<
" own:" << own <<
" nei:" << nei
348 "polyTopoChange::checkFace(const face&, const label"
349 ", const label, const label, const label)"
350 ) <<
"Cannot both have valid patchI and neighbour" <<
nl
352 <<
" faceI(-1 if added face):" << faceI
353 <<
" own:" << own <<
" nei:" << nei
361 "polyTopoChange::checkFace(const face&, const label"
362 ", const label, const label, const label)"
363 ) <<
"Owner cell label should be less than neighbour cell label"
366 <<
" faceI(-1 if added face):" << faceI
367 <<
" own:" << own <<
" nei:" << nei
372 if (f.size() < 3 ||
findIndex(f, -1) != -1)
376 "polyTopoChange::checkFace(const face&, const label"
377 ", const label, const label, const label)"
378 ) <<
"Illegal vertices in face"
381 <<
" faceI(-1 if added face):" << faceI
382 <<
" own:" << own <<
" nei:" << nei
385 if (faceI >= 0 && faceI < faces_.size() && faceRemoved(faceI))
389 "polyTopoChange::checkFace(const face&, const label"
390 ", const label, const label, const label)"
391 ) <<
"Face already marked for removal"
394 <<
" faceI(-1 if added face):" << faceI
395 <<
" own:" << own <<
" nei:" << nei
400 if (f[fp] < points_.size() && pointRemoved(f[fp]))
404 "polyTopoChange::checkFace(const face&, const label"
405 ", const label, const label, const label)"
406 ) <<
"Face uses removed vertices"
409 <<
" faceI(-1 if added face):" << faceI
410 <<
" own:" << own <<
" nei:" << nei
417 void Foam::polyTopoChange::makeCells
419 const label nActiveFaces,
424 cellFaces.setSize(2*nActiveFaces);
425 cellFaceOffsets.setSize(cellMap_.size() + 1);
432 for (label faceI = 0; faceI < nActiveFaces; faceI++)
434 nNbrs[faceOwner_[faceI]]++;
436 for (label faceI = 0; faceI < nActiveFaces; faceI++)
438 if (faceNeighbour_[faceI] >= 0)
440 nNbrs[faceNeighbour_[faceI]]++;
446 cellFaceOffsets[0] = 0;
449 cellFaceOffsets[cellI+1] = cellFaceOffsets[cellI] + nNbrs[cellI];
457 for (label faceI = 0; faceI < nActiveFaces; faceI++)
459 label cellI = faceOwner_[faceI];
461 cellFaces[cellFaceOffsets[cellI] + nNbrs[cellI]++] = faceI;
464 for (label faceI = 0; faceI < nActiveFaces; faceI++)
466 label cellI = faceNeighbour_[faceI];
470 cellFaces[cellFaceOffsets[cellI] + nNbrs[cellI]++] = faceI;
475 cellFaces.setSize(cellFaceOffsets[cellMap_.size()]);
481 void Foam::polyTopoChange::makeCellCells
483 const label nActiveFaces,
484 CompactListList<label>& cellCells
491 label nCellCells = 0;
495 for (label faceI = 0; faceI < nActiveFaces; faceI++)
497 if (faceNeighbour_[faceI] >= 0)
499 nNbrs[faceOwner_[faceI]]++;
500 nNbrs[faceNeighbour_[faceI]]++;
505 cellCells.
setSize(cellMap_.size(), nCellCells);
509 labelList& offsets = cellCells.offsets();
514 sumSize += nNbrs[cellI];
515 offsets[cellI] = sumSize;
523 for (label faceI = 0; faceI < nActiveFaces; faceI++)
525 label nei = faceNeighbour_[faceI];
529 label own = faceOwner_[faceI];
530 cellCells.m()[cellCells.index(own, nNbrs[own]++)] = nei;
531 cellCells.m()[cellCells.index(nei, nNbrs[nei]++)] = own;
539 Foam::label Foam::polyTopoChange::getCellOrder
541 const CompactListList<label>& cellCellAddressing,
545 const labelList& offsets = cellCellAddressing.offsets();
547 labelList newOrder(cellCellAddressing.size());
550 SLList<label> nextCell;
555 label cellInOrder = 0;
562 if (!cellRemoved(cellI) && visited.get(cellI) == 0)
565 nextCell.append(cellI);
574 label currentCell = nextCell.removeHead();
576 if (visited.get(currentCell) == 0)
578 visited.set(currentCell, 1);
581 newOrder[cellInOrder] = currentCell;
585 label i0 = (currentCell == 0 ? 0 : offsets[currentCell-1]);
586 label i1 = offsets[currentCell];
588 for (label i = i0; i < i1; i++)
590 label nbr = cellCellAddressing.m()[i];
592 if (!cellRemoved(nbr) && visited.get(nbr) == 0)
595 nextCell.append(nbr);
600 while (nextCell.size());
606 newOrder.setSize(cellInOrder);
609 oldToNew =
invert(cellCellAddressing.size(), newOrder);
618 void Foam::polyTopoChange::getFaceOrder
620 const label nActiveFaces,
629 oldToNew.setSize(faceOwner_.size());
637 label startOfCell = cellFaceOffsets[cellI];
638 label nFaces = cellFaceOffsets[cellI+1] - startOfCell;
641 SortableList<label> nbr(nFaces);
643 for (label i = 0; i < nFaces; i++)
645 label faceI = cellFaces[startOfCell + i];
647 label nbrCellI = faceNeighbour_[faceI];
649 if (faceI >= nActiveFaces)
654 else if (nbrCellI != -1)
657 if (nbrCellI == cellI)
659 nbrCellI = faceOwner_[faceI];
662 if (cellI < nbrCellI)
686 oldToNew[cellFaces[startOfCell + nbr.indices()[i]]] =
694 patchStarts.setSize(nPatches_);
696 patchSizes.setSize(nPatches_);
699 patchStarts[0] = newFaceI;
701 for (label faceI = 0; faceI < nActiveFaces; faceI++)
703 if (region_[faceI] >= 0)
705 patchSizes[region_[faceI]]++;
709 label faceI = patchStarts[0];
711 forAll(patchStarts, patchI)
713 patchStarts[patchI] = faceI;
714 faceI += patchSizes[patchI];
725 for (label faceI = 0; faceI < nActiveFaces; faceI++)
727 if (region_[faceI] >= 0)
729 oldToNew[faceI] = workPatchStarts[region_[faceI]]++;
734 for (label faceI = nActiveFaces; faceI < oldToNew.size(); faceI++)
736 oldToNew[faceI] = faceI;
742 if (oldToNew[faceI] == -1)
746 "polyTopoChange::getFaceOrder"
747 "(const label, const labelList&, const labelList&)"
749 ) <<
"Did not determine new position"
750 <<
" for face " << faceI
758 void Foam::polyTopoChange::reorderCompactFaces
765 faces_.setCapacity(newSize);
768 region_.setCapacity(newSize);
771 faceOwner_.setCapacity(newSize);
773 reorder(oldToNew, faceNeighbour_);
774 faceNeighbour_.setCapacity(newSize);
778 faceMap_.setCapacity(newSize);
780 renumberReverseMap(oldToNew, reverseFaceMap_);
782 renumberKey(oldToNew, faceFromPoint_);
783 renumberKey(oldToNew, faceFromEdge_);
785 flipFaceFlux_.setCapacity(newSize);
788 faceZone_.setCapacity(newSize);
791 faceZoneFlip_.setCapacity(newSize);
799 void Foam::polyTopoChange::compact
801 const bool orderCells,
802 const bool orderPoints,
803 label& nInternalPoints,
810 reversePointMap_.shrink();
816 faceNeighbour_.shrink();
818 reverseFaceMap_.shrink();
822 reverseCellMap_.shrink();
827 label nActivePoints = 0;
829 labelList localPointMap(points_.size(), -1);
834 nInternalPoints = -1;
838 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
840 localPointMap[pointI] = newPointI++;
843 nActivePoints = newPointI;
849 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
861 && faceOwner_[faceI] >= 0
862 && faceNeighbour_[faceI] < 0
866 const face& f = faces_[faceI];
870 label pointI = f[fp];
872 if (localPointMap[pointI] == -1)
877 || retiredPoints_.found(pointI)
881 <<
"Removed or retired point " << pointI
883 <<
" at position " << faceI <<
endl
884 <<
"Probably face has not been adapted for"
887 localPointMap[pointI] = newPointI++;
893 label nBoundaryPoints = newPointI;
894 nInternalPoints = nActivePoints - nBoundaryPoints;
897 forAll(localPointMap, pointI)
899 if (localPointMap[pointI] != -1)
901 localPointMap[pointI] += nInternalPoints;
913 && faceOwner_[faceI] >= 0
914 && faceNeighbour_[faceI] >= 0
918 const face& f = faces_[faceI];
922 label pointI = f[fp];
924 if (localPointMap[pointI] == -1)
929 || retiredPoints_.found(pointI)
933 <<
"Removed or retired point " << pointI
935 <<
" at position " << faceI <<
endl
936 <<
"Probably face has not been adapted for"
939 localPointMap[pointI] = newPointI++;
945 if (newPointI != nInternalPoints)
950 newPointI = nActivePoints;
955 localPointMap[iter.key()] = newPointI++;
961 Pout<<
"Points : active:" << nActivePoints
962 <<
" removed:" << points_.size()-newPointI <<
endl;
965 reorder(localPointMap, points_);
966 points_.setCapacity(newPointI);
969 reorder(localPointMap, pointMap_);
970 pointMap_.setCapacity(newPointI);
971 renumberReverseMap(localPointMap, reversePointMap_);
973 reorder(localPointMap, pointZone_);
974 pointZone_.setCapacity(newPointI);
975 renumber(localPointMap, retiredPoints_);
980 face& f = faces_[faceI];
983 renumberCompact(localPointMap, f);
985 if (!faceRemoved(faceI) && f.size() < 3)
988 <<
"Created illegal face " << f
990 <<
" at position:" << faceI
991 <<
" when filtering removed points"
1000 labelList localFaceMap(faces_.size(), -1);
1005 if (!faceRemoved(faceI) && faceOwner_[faceI] >= 0)
1007 localFaceMap[faceI] = newFaceI++;
1010 nActiveFaces_ = newFaceI;
1014 if (!faceRemoved(faceI) && faceOwner_[faceI] < 0)
1017 localFaceMap[faceI] = newFaceI++;
1023 Pout<<
"Faces : active:" << nActiveFaces_
1024 <<
" removed:" << faces_.size()-newFaceI <<
endl;
1028 reorderCompactFaces(newFaceI, localFaceMap);
1039 CompactListList<label> cellCells;
1040 makeCellCells(nActiveFaces_, cellCells);
1043 newCellI = getCellOrder(cellCells, localCellMap);
1048 localCellMap.setSize(cellMap_.size());
1054 if (!cellRemoved(cellI))
1056 localCellMap[cellI] = newCellI++;
1063 Pout<<
"Cells : active:" << newCellI
1064 <<
" removed:" << cellMap_.size()-newCellI <<
endl;
1068 if (orderCells || (newCellI != cellMap_.size()))
1070 reorder(localCellMap, cellMap_);
1071 cellMap_.setCapacity(newCellI);
1072 renumberReverseMap(localCellMap, reverseCellMap_);
1074 reorder(localCellMap, cellZone_);
1075 cellZone_.setCapacity(newCellI);
1077 renumberKey(localCellMap, cellFromPoint_);
1078 renumberKey(localCellMap, cellFromEdge_);
1079 renumberKey(localCellMap, cellFromFace_);
1083 forAll(faceOwner_, faceI)
1085 label own = faceOwner_[faceI];
1086 label nei = faceNeighbour_[faceI];
1091 faceOwner_[faceI] = localCellMap[own];
1096 faceNeighbour_[faceI] = localCellMap[nei];
1101 faceNeighbour_[faceI] >= 0
1102 && faceNeighbour_[faceI] < faceOwner_[faceI]
1105 faces_[faceI] = faces_[faceI].reverseFace();
1106 Swap(faceOwner_[faceI], faceNeighbour_[faceI]);
1107 flipFaceFlux_[faceI] =
1109 flipFaceFlux_[faceI]
1113 faceZoneFlip_[faceI] =
1115 faceZoneFlip_[faceI]
1125 faceNeighbour_[faceI] = localCellMap[nei];
1136 makeCells(nActiveFaces_, cellFaces, cellFaceOffsets);
1152 reorderCompactFaces(localFaceMap.size(), localFaceMap);
1166 const primitiveMesh& mesh,
1168 const bool internalFacesOnly
1175 label faceI = faceLabels[i];
1177 if (internalFacesOnly == mesh.isInternalFace(faceI))
1189 collectedFaces = faceLabels;
1193 collectedFaces.
setSize(nFaces);
1199 label faceI = faceLabels[i];
1201 if (internalFacesOnly == mesh.isInternalFace(faceI))
1203 collectedFaces[nFaces++] = faceI;
1208 return collectedFaces;
1214 void Foam::polyTopoChange::calcPatchPointMap
1216 const List<Map<label> >& oldPatchMeshPointMaps,
1221 patchPointMap.setSize(boundary.size());
1225 const labelList& meshPoints = boundary[patchI].meshPoints();
1227 const Map<label>& oldMeshPointMap = oldPatchMeshPointMaps[patchI];
1229 labelList& curPatchPointRnb = patchPointMap[patchI];
1231 curPatchPointRnb.
setSize(meshPoints.size());
1235 if (meshPoints[i] < pointMap_.size())
1240 pointMap_[meshPoints[i]]
1243 if (ozmpmIter != oldMeshPointMap.end())
1245 curPatchPointRnb[i] = ozmpmIter();
1249 curPatchPointRnb[i] = -1;
1254 curPatchPointRnb[i] = -1;
1261 void Foam::polyTopoChange::calcFaceInflationMaps
1263 const polyMesh& mesh,
1264 List<objectMap>& facesFromPoints,
1265 List<objectMap>& facesFromEdges,
1266 List<objectMap>& facesFromFaces
1272 facesFromPoints.setSize(faceFromPoint_.size());
1274 if (faceFromPoint_.size())
1276 label nFacesFromPoints = 0;
1281 label newFaceI = iter.key();
1283 if (region_[newFaceI] == -1)
1286 facesFromPoints[nFacesFromPoints++] = objectMap
1292 mesh.pointFaces()[iter()],
1300 facesFromPoints[nFacesFromPoints++] = objectMap
1306 mesh.pointFaces()[iter()],
1318 facesFromEdges.setSize(faceFromEdge_.size());
1320 if (faceFromEdge_.size())
1322 label nFacesFromEdges = 0;
1327 label newFaceI = iter.key();
1329 if (region_[newFaceI] == -1)
1332 facesFromEdges[nFacesFromEdges++] = objectMap
1338 mesh.edgeFaces(iter()),
1346 facesFromEdges[nFacesFromEdges++] = objectMap
1352 mesh.edgeFaces(iter()),
1373 void Foam::polyTopoChange::calcCellInflationMaps
1375 const polyMesh& mesh,
1376 List<objectMap>& cellsFromPoints,
1377 List<objectMap>& cellsFromEdges,
1378 List<objectMap>& cellsFromFaces,
1379 List<objectMap>& cellsFromCells
1382 cellsFromPoints.setSize(cellFromPoint_.size());
1384 if (cellFromPoint_.size())
1386 label nCellsFromPoints = 0;
1391 cellsFromPoints[nCellsFromPoints++] = objectMap
1394 mesh.pointCells()[iter()]
1400 cellsFromEdges.setSize(cellFromEdge_.size());
1402 if (cellFromEdge_.size())
1404 label nCellsFromEdges = 0;
1409 cellsFromEdges[nCellsFromEdges++] = objectMap
1412 mesh.edgeCells()[iter()]
1418 cellsFromFaces.setSize(cellFromFace_.size());
1420 if (cellFromFace_.size())
1422 label nCellsFromFaces = 0;
1429 label oldFaceI = iter();
1431 if (mesh.isInternalFace(oldFaceI))
1433 twoCells[0] = mesh.faceOwner()[oldFaceI];
1434 twoCells[1] = mesh.faceNeighbour()[oldFaceI];
1435 cellsFromFaces[nCellsFromFaces++] = objectMap
1443 cellsFromFaces[nCellsFromFaces++] = objectMap
1446 labelList(1, mesh.faceOwner()[oldFaceI])
1465 void Foam::polyTopoChange::resetZones
1467 const polyMesh& mesh,
1477 pointZoneMap.setSize(mesh.pointZones().size());
1485 forAll(pointZone_, pointI)
1487 label zoneI = pointZone_[pointI];
1489 if (zoneI >= pointZones.size())
1493 "resetZones(const polyMesh&, polyMesh&, labelListList&"
1494 "labelListList&, labelListList&)"
1495 ) <<
"Illegal zoneID " << zoneI <<
" for point "
1496 << pointI <<
" coord " << mesh.points()[pointI]
1509 forAll(addressing, zoneI)
1511 addressing[zoneI].setSize(
nPoints[zoneI]);
1515 forAll(pointZone_, pointI)
1517 label zoneI = pointZone_[pointI];
1520 addressing[zoneI][
nPoints[zoneI]++] = pointI;
1524 forAll(addressing, zoneI)
1531 forAll(addressing, zoneI)
1533 const pointZone& oldZone = pointZones[zoneI];
1534 const labelList& newZoneAddr = addressing[zoneI];
1536 labelList& curPzRnb = pointZoneMap[zoneI];
1537 curPzRnb.
setSize(newZoneAddr.size());
1541 if (newZoneAddr[i] < pointMap_.size())
1543 curPzRnb[i] = oldZone.whichPoint(pointMap_[newZoneAddr[i]]);
1553 newMesh.pointZones().clearAddressing();
1554 forAll(newMesh.pointZones(), zoneI)
1558 Pout<<
"pointZone:" << zoneI
1559 <<
" name:" << newMesh.pointZones()[zoneI].
name()
1560 <<
" size:" << addressing[zoneI].size()
1564 newMesh.pointZones()[zoneI] = addressing[zoneI];
1572 faceZoneFaceMap.setSize(mesh.faceZones().size());
1580 label zoneI = faceZone_[faceI];
1582 if (zoneI >= faceZones.size())
1586 "resetZones(const polyMesh&, polyMesh&, labelListList&"
1587 "labelListList&, labelListList&)"
1588 ) <<
"Illegal zoneID " << zoneI <<
" for face "
1601 forAll(addressing, zoneI)
1603 addressing[zoneI].setSize(nFaces[zoneI]);
1604 flipMode[zoneI].setSize(nFaces[zoneI]);
1610 label zoneI = faceZone_[faceI];
1614 label index = nFaces[zoneI]++;
1615 addressing[zoneI][index] = faceI;
1616 flipMode[zoneI][index] = faceZoneFlip_[faceI];
1620 forAll(addressing, zoneI)
1625 labelList newAddressing(addressing[zoneI].size());
1628 newAddressing[i] = addressing[zoneI][newToOld[i]];
1630 addressing[zoneI].transfer(newAddressing);
1633 boolList newFlipMode(flipMode[zoneI].size());
1636 newFlipMode[i] = flipMode[zoneI][newToOld[i]];
1638 flipMode[zoneI].transfer(newFlipMode);
1644 forAll(addressing, zoneI)
1646 const faceZone& oldZone = faceZones[zoneI];
1647 const labelList& newZoneAddr = addressing[zoneI];
1649 labelList& curFzFaceRnb = faceZoneFaceMap[zoneI];
1651 curFzFaceRnb.
setSize(newZoneAddr.size());
1655 if (newZoneAddr[i] < faceMap_.size())
1658 oldZone.whichFace(faceMap_[newZoneAddr[i]]);
1662 curFzFaceRnb[i] = -1;
1669 newMesh.faceZones().clearAddressing();
1670 forAll(newMesh.faceZones(), zoneI)
1674 Pout<<
"faceZone:" << zoneI
1675 <<
" name:" << newMesh.faceZones()[zoneI].
name()
1676 <<
" size:" << addressing[zoneI].size()
1680 newMesh.faceZones()[zoneI].resetAddressing
1692 cellZoneMap.setSize(mesh.cellZones().size());
1700 label zoneI = cellZone_[cellI];
1702 if (zoneI >= cellZones.size())
1706 "resetZones(const polyMesh&, polyMesh&, labelListList&"
1707 "labelListList&, labelListList&)"
1708 ) <<
"Illegal zoneID " << zoneI <<
" for cell "
1719 forAll(addressing, zoneI)
1721 addressing[zoneI].setSize(nCells[zoneI]);
1727 label zoneI = cellZone_[cellI];
1731 addressing[zoneI][nCells[zoneI]++] = cellI;
1735 forAll(addressing, zoneI)
1742 forAll(addressing, zoneI)
1744 const cellZone& oldZone = cellZones[zoneI];
1745 const labelList& newZoneAddr = addressing[zoneI];
1747 labelList& curCellRnb = cellZoneMap[zoneI];
1749 curCellRnb.
setSize(newZoneAddr.size());
1753 if (newZoneAddr[i] < cellMap_.size())
1756 oldZone.whichCell(cellMap_[newZoneAddr[i]]);
1766 newMesh.cellZones().clearAddressing();
1767 forAll(newMesh.cellZones(), zoneI)
1771 Pout<<
"cellZone:" << zoneI
1772 <<
" name:" << newMesh.cellZones()[zoneI].
name()
1773 <<
" size:" << addressing[zoneI].size()
1777 newMesh.cellZones()[zoneI] = addressing[zoneI];
1783 void Foam::polyTopoChange::calcFaceZonePointMap
1785 const polyMesh& mesh,
1786 const List<Map<label> >& oldFaceZoneMeshPointMaps,
1792 faceZonePointMap.
setSize(faceZones.size());
1796 const faceZone& newZone = faceZones[zoneI];
1798 const labelList& newZoneMeshPoints = newZone().meshPoints();
1800 const Map<label>& oldZoneMeshPointMap = oldFaceZoneMeshPointMaps[zoneI];
1802 labelList& curFzPointRnb = faceZonePointMap[zoneI];
1804 curFzPointRnb.
setSize(newZoneMeshPoints.size());
1806 forAll(newZoneMeshPoints, pointI)
1808 if (newZoneMeshPoints[pointI] < pointMap_.size())
1811 oldZoneMeshPointMap.
find
1813 pointMap_[newZoneMeshPoints[pointI]]
1816 if (ozmpmIter != oldZoneMeshPointMap.end())
1818 curFzPointRnb[pointI] = ozmpmIter();
1822 curFzPointRnb[pointI] = -1;
1827 curFzPointRnb[pointI] = -1;
1840 face newF(f.size());
1844 label fp1 = (fp + nPos) % f.size();
1858 void Foam::polyTopoChange::reorderCoupledFaces
1860 const bool syncParallel,
1861 const polyBoundaryMesh& boundary,
1877 if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1879 boundary[patchI].initOrder
1897 bool anyChanged =
false;
1901 if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1903 labelList patchFaceMap(patchSizes[patchI], -1);
1904 labelList patchFaceRotation(patchSizes[patchI], 0);
1906 bool changed = boundary[patchI].order
1925 label start = patchStarts[patchI];
1927 forAll(patchFaceMap, patchFaceI)
1929 oldToNew[patchFaceI + start] =
1930 start + patchFaceMap[patchFaceI];
1933 forAll(patchFaceRotation, patchFaceI)
1935 rotation[patchFaceI + start] =
1936 patchFaceRotation[patchFaceI];
1946 reduce(anyChanged, orOp<bool>());
1952 reorderCompactFaces(oldToNew.size(), oldToNew);
1957 if (rotation[faceI] != 0)
1959 faces_[faceI] = rotateFace(faces_[faceI], rotation[faceI]);
1966 void Foam::polyTopoChange::compactAndReorder
1968 const polyMesh& mesh,
1969 const bool syncParallel,
1970 const bool orderCells,
1971 const bool orderPoints,
1973 label& nInternalPoints,
1977 List<objectMap>& pointsFromPoints,
1978 List<objectMap>& facesFromPoints,
1979 List<objectMap>& facesFromEdges,
1980 List<objectMap>& facesFromFaces,
1981 List<objectMap>& cellsFromPoints,
1982 List<objectMap>& cellsFromEdges,
1983 List<objectMap>& cellsFromFaces,
1984 List<objectMap>& cellsFromCells,
1985 List<Map<label> >& oldPatchMeshPointMaps,
1988 List<Map<label> >& oldFaceZoneMeshPointMaps
1991 if (mesh.boundaryMesh().size() != nPatches_)
1994 <<
"polyTopoChange was constructed with a mesh with "
1995 << nPatches_ <<
" patches." <<
endl
1996 <<
"The mesh now provided has a different number of patches "
1997 << mesh.boundaryMesh().size()
1998 <<
" which is illegal" <<
endl
2004 compact(orderCells, orderPoints, nInternalPoints, patchSizes, patchStarts);
2008 newPoints.transfer(points_);
2014 mesh.boundaryMesh(),
2040 calcFaceInflationMaps
2048 calcCellInflationMaps
2059 faceFromPoint_.clearStorage();
2060 faceFromEdge_.clearStorage();
2062 cellFromPoint_.clearStorage();
2063 cellFromEdge_.clearStorage();
2064 cellFromFace_.clearStorage();
2068 const polyBoundaryMesh& boundary = mesh.boundaryMesh();
2071 oldPatchMeshPointMaps.setSize(boundary.size());
2072 oldPatchNMeshPoints.setSize(boundary.size());
2073 oldPatchStarts.setSize(boundary.size());
2078 oldPatchMeshPointMaps[patchI] = boundary[patchI].meshPointMap();
2079 oldPatchNMeshPoints[patchI] = boundary[patchI].meshPoints().size();
2080 oldPatchStarts[patchI] = boundary[patchI].start();
2086 oldFaceZoneMeshPointMaps.setSize(mesh.faceZones().size());
2088 forAll(mesh.faceZones(), zoneI)
2090 const faceZone& oldZone = mesh.faceZones()[zoneI];
2092 oldFaceZoneMeshPointMaps[zoneI] = oldZone().meshPointMap();
2103 nPatches_(nPatches),
2106 reversePointMap_(0),
2141 reversePointMap_(0),
2178 points_.clearStorage();
2179 pointMap_.clearStorage();
2180 reversePointMap_.clearStorage();
2181 pointZone_.clearStorage();
2182 retiredPoints_.clearStorage();
2184 faces_.clearStorage();
2185 region_.clearStorage();
2186 faceOwner_.clearStorage();
2187 faceNeighbour_.clearStorage();
2188 faceMap_.clearStorage();
2189 reverseFaceMap_.clearStorage();
2190 faceFromPoint_.clearStorage();
2191 faceFromEdge_.clearStorage();
2192 flipFaceFlux_.clearStorage();
2193 faceZone_.clearStorage();
2194 faceZoneFlip_.clearStorage();
2197 cellMap_.clearStorage();
2198 reverseCellMap_.clearStorage();
2199 cellZone_.clearStorage();
2200 cellFromPoint_.clearStorage();
2201 cellFromEdge_.clearStorage();
2202 cellFromFace_.clearStorage();
2215 label maxRegion = nPatches_ - 1;
2218 maxRegion =
max(maxRegion, patchMap[i]);
2220 nPatches_ = maxRegion + 1;
2229 points_.setCapacity(points_.size() + points.
size());
2230 pointMap_.setCapacity(pointMap_.size() + points.
size());
2231 reversePointMap_.setCapacity(reversePointMap_.size() + points.
size());
2232 pointZone_.setCapacity(pointZone_.size() + points.
size());
2237 forAll(pointZones, zoneI)
2243 newZoneID[pointLabels[j]] = pointZoneMap[zoneI];
2248 for (label pointI = 0; pointI < mesh.
nPoints(); pointI++)
2268 label nAllCells = mesh.
nCells();
2270 cellMap_.setCapacity(cellMap_.size() + nAllCells);
2271 reverseCellMap_.setCapacity(reverseCellMap_.size() + nAllCells);
2272 cellFromPoint_.resize(cellFromPoint_.size() + nAllCells/100);
2273 cellFromEdge_.resize(cellFromEdge_.size() + nAllCells/100);
2274 cellFromFace_.resize(cellFromFace_.size() + nAllCells/100);
2275 cellZone_.setCapacity(cellZone_.size() + nAllCells);
2283 const labelList& cellLabels = cellZones[zoneI];
2287 label cellI = cellLabels[j];
2289 if (newZoneID[cellI] != -1)
2293 "polyTopoChange::addMesh"
2294 "(const polyMesh&, const labelList&,"
2295 "const labelList&, const labelList&,"
2297 ) <<
"Cell:" << cellI
2299 <<
" is in two zones:"
2300 << cellZones[newZoneID[cellI]].
name()
2301 <<
" and " << cellZones[zoneI].
name() <<
endl
2302 <<
" This is not supported."
2303 <<
" Continuing with first zone only." <<
endl;
2307 newZoneID[cellI] = cellZoneMap[zoneI];
2313 for (label cellI = 0; cellI < nAllCells; cellI++)
2316 addCell(-1, -1, -1, cellI, newZoneID[cellI]);
2331 faces_.setCapacity(faces_.size() + nAllFaces);
2332 region_.setCapacity(region_.size() + nAllFaces);
2333 faceOwner_.setCapacity(faceOwner_.size() + nAllFaces);
2334 faceNeighbour_.setCapacity(faceNeighbour_.size() + nAllFaces);
2335 faceMap_.setCapacity(faceMap_.size() + nAllFaces);
2336 reverseFaceMap_.setCapacity(reverseFaceMap_.size() + nAllFaces);
2337 faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/100);
2338 faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/100);
2339 flipFaceFlux_.setCapacity(flipFaceFlux_.size() + nAllFaces);
2340 faceZone_.setCapacity(faceZone_.size() + nAllFaces);
2341 faceZoneFlip_.setCapacity(faceZoneFlip_.size() + nAllFaces);
2346 boolList zoneFlip(nAllFaces,
false);
2350 const labelList& faceLabels = faceZones[zoneI];
2351 const boolList& flipMap = faceZones[zoneI].flipMap();
2355 newZoneID[faceLabels[j]] = faceZoneMap[zoneI];
2356 zoneFlip[faceLabels[j]] = flipMap[j];
2369 faceNeighbour[faceI],
2385 if (pp.
start() != faces_.size())
2389 "polyTopoChange::polyTopoChange"
2390 "(const polyMesh& mesh, const bool strict)"
2392 <<
"Patch " << pp.
name() <<
" starts at " << pp.
start()
2394 <<
"Current face counter at " << faces_.size() <<
endl
2395 <<
"Are patches in incremental order?"
2400 label faceI = pp.
start() + patchFaceI;
2428 points_.setCapacity(nPoints);
2429 pointMap_.setCapacity(nPoints);
2430 reversePointMap_.setCapacity(nPoints);
2431 pointZone_.setCapacity(nPoints);
2433 faces_.setCapacity(nFaces);
2434 region_.setCapacity(nFaces);
2435 faceOwner_.setCapacity(nFaces);
2436 faceNeighbour_.setCapacity(nFaces);
2437 faceMap_.setCapacity(nFaces);
2438 reverseFaceMap_.setCapacity(nFaces);
2439 faceFromPoint_.resize(faceFromPoint_.size() + nFaces/100);
2440 faceFromEdge_.resize(faceFromEdge_.size() + nFaces/100);
2441 flipFaceFlux_.setCapacity(nFaces);
2442 faceZone_.setCapacity(nFaces);
2443 faceZoneFlip_.setCapacity(nFaces);
2445 cellMap_.setCapacity(nCells);
2446 reverseCellMap_.setCapacity(nCells);
2447 cellFromPoint_.resize(cellFromPoint_.size() + nCells/100);
2448 cellFromEdge_.resize(cellFromEdge_.size() + nCells/100);
2449 cellFromFace_.resize(cellFromFace_.size() + nCells/100);
2450 cellZone_.setCapacity(nCells);
2456 if (isType<polyAddPoint>(action))
2458 const polyAddPoint& pap = refCast<const polyAddPoint>(action);
2468 else if (isType<polyModifyPoint>(action))
2482 else if (isType<polyRemovePoint>(action))
2490 else if (isType<polyAddFace>(action))
2492 const polyAddFace& paf = refCast<const polyAddFace>(action);
2508 else if (isType<polyModifyFace>(action))
2510 const polyModifyFace& pmf = refCast<const polyModifyFace>(action);
2526 else if (isType<polyRemoveFace>(action))
2528 const polyRemoveFace& prf = refCast<const polyRemoveFace>(action);
2534 else if (isType<polyAddCell>(action))
2536 const polyAddCell& pac = refCast<const polyAddCell>(action);
2547 else if (isType<polyModifyCell>(action))
2549 const polyModifyCell& pmc = refCast<const polyModifyCell>(action);
2553 modifyCell(pmc.
cellID(), -1);
2562 else if (isType<polyRemoveCell>(action))
2564 const polyRemoveCell& prc = refCast<const polyRemoveCell>(action);
2574 "label polyTopoChange::setAction(const topoAction& action)"
2575 ) <<
"Unknown type of topoChange: " << action.type()
2587 const label masterPointID,
2592 label pointI = points_.size();
2595 pointMap_.append(masterPointID);
2596 reversePointMap_.append(pointI);
2597 pointZone_.append(zoneID);
2601 retiredPoints_.insert(pointI);
2612 const label newZoneID,
2616 if (pointI < 0 || pointI >= points_.size())
2620 "polyTopoChange::modifyPoint(const label, const point&)"
2621 ) <<
"illegal point label " << pointI <<
endl
2622 <<
"Valid point labels are 0 .. " << points_.size()-1
2625 if (pointRemoved(pointI) || pointMap_[pointI] == -1)
2629 "polyTopoChange::modifyPoint(const label, const point&)"
2630 ) <<
"point " << pointI <<
" already marked for removal"
2633 points_[pointI] = pt;
2634 pointZone_[pointI] = newZoneID;
2638 retiredPoints_.erase(pointI);
2642 retiredPoints_.insert(pointI);
2649 if (newPoints.
size() != points_.size())
2651 FatalErrorIn(
"polyTopoChange::movePoints(const pointField&)")
2652 <<
"illegal pointField size." <<
endl
2653 <<
"Size:" << newPoints.
size() <<
endl
2654 <<
"Points in mesh:" << points_.size()
2660 points_[pointI] = newPoints[pointI];
2668 const label mergePointI
2671 if (pointI < 0 || pointI >= points_.size())
2673 FatalErrorIn(
"polyTopoChange::removePoint(const label, const label)")
2674 <<
"illegal point label " << pointI <<
endl
2675 <<
"Valid point labels are 0 .. " << points_.size()-1
2682 && (pointRemoved(pointI) || pointMap_[pointI] == -1)
2685 FatalErrorIn(
"polyTopoChange::removePoint(const label, const label)")
2686 <<
"point " << pointI <<
" already marked for removal" <<
nl
2687 <<
"Point:" << points_[pointI] <<
" pointMap:" << pointMap_[pointI]
2691 if (pointI == mergePointI)
2693 FatalErrorIn(
"polyTopoChange::removePoint(const label, const label)")
2694 <<
"Cannot remove/merge point " << pointI <<
" onto itself."
2698 points_[pointI] = greatPoint;
2699 pointMap_[pointI] = -1;
2700 if (mergePointI >= 0)
2702 reversePointMap_[pointI] = -mergePointI-2;
2706 reversePointMap_[pointI] = -1;
2708 pointZone_[pointI] = -1;
2709 retiredPoints_.erase(pointI);
2718 const label masterPointID,
2719 const label masterEdgeID,
2720 const label masterFaceID,
2721 const bool flipFaceFlux,
2722 const label patchID,
2730 checkFace(f, -1, own, nei, patchID, zoneID);
2733 label faceI = faces_.size();
2736 region_.append(patchID);
2737 faceOwner_.append(own);
2738 faceNeighbour_.append(nei);
2740 if (masterPointID >= 0)
2742 faceMap_.append(-1);
2743 faceFromPoint_.insert(faceI, masterPointID);
2745 else if (masterEdgeID >= 0)
2747 faceMap_.append(-1);
2748 faceFromEdge_.insert(faceI, masterEdgeID);
2750 else if (masterFaceID >= 0)
2752 faceMap_.append(masterFaceID);
2761 faceMap_.append(-1);
2763 reverseFaceMap_.append(faceI);
2765 flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2766 faceZone_.append(zoneID);
2767 faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2779 const bool flipFaceFlux,
2780 const label patchID,
2788 checkFace(f, faceI, own, nei, patchID, zoneID);
2792 faceOwner_[faceI] = own;
2793 faceNeighbour_[faceI] = nei;
2794 region_[faceI] = patchID;
2796 flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2797 faceZone_[faceI] = zoneID;
2798 faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2804 if (faceI < 0 || faceI >= faces_.size())
2806 FatalErrorIn(
"polyTopoChange::removeFace(const label, const label)")
2807 <<
"illegal face label " << faceI <<
endl
2808 <<
"Valid face labels are 0 .. " << faces_.size()-1
2815 && (faceRemoved(faceI) || faceMap_[faceI] == -1)
2818 FatalErrorIn(
"polyTopoChange::removeFace(const label, const label)")
2820 <<
" already marked for removal"
2824 faces_[faceI].setSize(0);
2825 region_[faceI] = -1;
2826 faceOwner_[faceI] = -1;
2827 faceNeighbour_[faceI] = -1;
2828 faceMap_[faceI] = -1;
2829 if (mergeFaceI >= 0)
2831 reverseFaceMap_[faceI] = -mergeFaceI-2;
2835 reverseFaceMap_[faceI] = -1;
2837 faceFromEdge_.erase(faceI);
2838 faceFromPoint_.erase(faceI);
2839 flipFaceFlux_[faceI] = 0;
2840 faceZone_[faceI] = -1;
2841 faceZoneFlip_[faceI] = 0;
2847 const label masterPointID,
2848 const label masterEdgeID,
2849 const label masterFaceID,
2850 const label masterCellID,
2854 label cellI = cellMap_.size();
2856 if (masterPointID >= 0)
2858 cellMap_.append(-1);
2859 cellFromPoint_.insert(cellI, masterPointID);
2861 else if (masterEdgeID >= 0)
2863 cellMap_.append(-1);
2864 cellFromEdge_.insert(cellI, masterEdgeID);
2866 else if (masterFaceID >= 0)
2868 cellMap_.append(-1);
2869 cellFromFace_.insert(cellI, masterFaceID);
2873 cellMap_.append(masterCellID);
2875 reverseCellMap_.append(cellI);
2876 cellZone_.append(zoneID);
2888 cellZone_[cellI] = zoneID;
2894 if (cellI < 0 || cellI >= cellMap_.size())
2896 FatalErrorIn(
"polyTopoChange::removeCell(const label, const label)")
2897 <<
"illegal cell label " << cellI <<
endl
2898 <<
"Valid cell labels are 0 .. " << cellMap_.size()-1
2902 if (strict_ && cellMap_[cellI] == -2)
2904 FatalErrorIn(
"polyTopoChange::removeCell(const label, const label)")
2906 <<
" already marked for removal"
2910 cellMap_[cellI] = -2;
2911 if (mergeCellI >= 0)
2913 reverseCellMap_[cellI] = -mergeCellI-2;
2917 reverseCellMap_[cellI] = -1;
2919 cellFromPoint_.erase(cellI);
2920 cellFromEdge_.erase(cellI);
2921 cellFromFace_.erase(cellI);
2922 cellZone_[cellI] = -1;
2930 const bool syncParallel,
2931 const bool orderCells,
2932 const bool orderPoints
2937 Pout<<
"polyTopoChange::changeMesh"
2938 <<
"(polyMesh&, const bool, const bool, const bool, const bool)"
2944 Pout<<
"Old mesh:" <<
nl;
2945 writeMeshStats(mesh,
Pout);
2951 label nInternalPoints;
2990 oldPatchMeshPointMaps,
2991 oldPatchNMeshPoints,
2993 oldFaceZoneMeshPointMaps
2996 const label nOldPoints(mesh.
nPoints());
2997 const label nOldFaces(mesh.
nFaces());
2998 const label nOldCells(mesh.
nCells());
3012 forAll(pointMap_, newPointI)
3014 label oldPointI = pointMap_[newPointI];
3018 renumberedMeshPoints[newPointI] = mesh.
points()[oldPointI];
3031 faceNeighbour_.xfer(),
3047 faceNeighbour_.xfer(),
3057 retiredPoints_.clearStorage();
3058 region_.clearStorage();
3065 label nAdd, nInflate, nMerge, nRemove;
3066 countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3068 <<
" added(from point):" << nAdd
3069 <<
" added(from nothing):" << nInflate
3070 <<
" merged(into other point):" << nMerge
3071 <<
" removed:" << nRemove
3074 countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3076 <<
" added(from face):" << nAdd
3077 <<
" added(inflated):" << nInflate
3078 <<
" merged(into other face):" << nMerge
3079 <<
" removed:" << nRemove
3082 countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3084 <<
" added(from cell):" << nAdd
3085 <<
" added(inflated):" << nInflate
3086 <<
" merged(into other cell):" << nMerge
3087 <<
" removed:" << nRemove
3094 Pout<<
"New mesh:" <<
nl;
3095 writeMeshStats(mesh,
Pout);
3109 resetZones(mesh, mesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3113 pointZone_.clearStorage();
3114 faceZone_.clearStorage();
3115 faceZoneFlip_.clearStorage();
3116 cellZone_.clearStorage();
3126 oldPatchMeshPointMaps,
3133 calcFaceZonePointMap(mesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3135 labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3176 oldPatchNMeshPoints,
3191 const bool syncParallel,
3192 const bool orderCells,
3193 const bool orderPoints
3198 Pout<<
"polyTopoChange::changeMesh"
3199 <<
"(autoPtr<fvMesh>&, const IOobject&, const fvMesh&"
3200 <<
", const bool, const bool, const bool)"
3206 Pout<<
"Old mesh:" <<
nl;
3207 writeMeshStats(mesh,
Pout);
3213 label nInternalPoints;
3253 oldPatchMeshPointMaps,
3254 oldPatchNMeshPoints,
3256 oldFaceZoneMeshPointMaps
3259 const label nOldPoints(mesh.
nPoints());
3260 const label nOldFaces(mesh.
nFaces());
3261 const label nOldCells(mesh.
nCells());
3275 faceNeighbour_.xfer()
3278 fvMesh& newMesh = newMeshPtr();
3282 retiredPoints_.clearStorage();
3283 region_.clearStorage();
3290 label nAdd, nInflate, nMerge, nRemove;
3291 countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3293 <<
" added(from point):" << nAdd
3294 <<
" added(from nothing):" << nInflate
3295 <<
" merged(into other point):" << nMerge
3296 <<
" removed:" << nRemove
3299 countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3301 <<
" added(from face):" << nAdd
3302 <<
" added(inflated):" << nInflate
3303 <<
" merged(into other face):" << nMerge
3304 <<
" removed:" << nRemove
3307 countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3309 <<
" added(from cell):" << nAdd
3310 <<
" added(inflated):" << nInflate
3311 <<
" merged(into other cell):" << nMerge
3312 <<
" removed:" << nRemove
3323 forAll(oldPatches, patchI)
3325 newBoundary[patchI] = oldPatches[patchI].
clone
3348 oldPointZones[i].
name(),
3363 oldFaceZones[i].
name(),
3379 oldCellZones[i].
name(),
3387 newMesh.
addZones(pZonePtrs, fZonePtrs, cZonePtrs);
3396 resetZones(mesh, newMesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3400 pointZone_.clearStorage();
3401 faceZone_.clearStorage();
3402 faceZoneFlip_.clearStorage();
3403 cellZone_.clearStorage();
3412 oldPatchMeshPointMaps,
3419 calcFaceZonePointMap(newMesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3423 Pout<<
"New mesh:" <<
nl;
3424 writeMeshStats(mesh,
Pout);
3427 labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3468 oldPatchNMeshPoints,