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);
787 renumberKey(oldToNew, faceZone_);
788 inplaceReorder(oldToNew, faceZoneFlip_);
789 faceZoneFlip_.setCapacity(newSize);
793 // Compact all and orders points and faces:
794 // - points into internal followed by external points
795 // - internalfaces upper-triangular
796 // - externalfaces after internal ones.
797 void Foam::polyTopoChange::compact
799 const bool orderCells,
800 const bool orderPoints,
801 label& nInternalPoints,
802 labelList& patchSizes,
803 labelList& patchStarts
808 reversePointMap_.shrink();
813 faceNeighbour_.shrink();
815 reverseFaceMap_.shrink();
818 reverseCellMap_.shrink();
823 label nActivePoints = 0;
825 labelList localPointMap(points_.size(), -1);
830 nInternalPoints = -1;
832 forAll(points_, pointI)
834 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
836 localPointMap[pointI] = newPointI++;
839 nActivePoints = newPointI;
843 forAll(points_, pointI)
845 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
851 // Mark boundary points
852 forAll(faceOwner_, faceI)
857 && faceOwner_[faceI] >= 0
858 && faceNeighbour_[faceI] < 0
861 // Valid boundary face
862 const face& f = faces_[faceI];
866 label pointI = f[fp];
868 if (localPointMap[pointI] == -1)
873 || retiredPoints_.found(pointI)
876 FatalErrorIn("polyTopoChange::compact(..)")
877 << "Removed or retired point " << pointI
879 << " at position " << faceI << endl
880 << "Probably face has not been adapted for"
881 << " removed points." << abort(FatalError);
883 localPointMap[pointI] = newPointI++;
889 label nBoundaryPoints = newPointI;
890 nInternalPoints = nActivePoints - nBoundaryPoints;
892 // Move the boundary addressing up
893 forAll(localPointMap, pointI)
895 if (localPointMap[pointI] != -1)
897 localPointMap[pointI] += nInternalPoints;
903 // Mark internal points
904 forAll(faceOwner_, faceI)
909 && faceOwner_[faceI] >= 0
910 && faceNeighbour_[faceI] >= 0
913 // Valid internal face
914 const face& f = faces_[faceI];
918 label pointI = f[fp];
920 if (localPointMap[pointI] == -1)
925 || retiredPoints_.found(pointI)
928 FatalErrorIn("polyTopoChange::compact(..)")
929 << "Removed or retired point " << pointI
931 << " at position " << faceI << endl
932 << "Probably face has not been adapted for"
933 << " removed points." << abort(FatalError);
935 localPointMap[pointI] = newPointI++;
941 if (newPointI != nInternalPoints)
943 FatalErrorIn("polyTopoChange::compact(..)")
944 << "Problem." << abort(FatalError);
946 newPointI = nActivePoints;
949 forAllConstIter(labelHashSet, retiredPoints_, iter)
951 localPointMap[iter.key()] = newPointI++;
957 Pout<< "Points : active:" << nActivePoints
958 << " removed:" << points_.size()-newPointI << endl;
961 reorder(localPointMap, points_);
962 points_.setCapacity(newPointI);
965 reorder(localPointMap, pointMap_);
966 pointMap_.setCapacity(newPointI);
967 renumberReverseMap(localPointMap, reversePointMap_);
969 renumberKey(localPointMap, pointZone_);
970 renumber(localPointMap, retiredPoints_);
972 // Use map to relabel face vertices
973 forAll(faces_, faceI)
975 face& f = faces_[faceI];
978 renumberCompact(localPointMap, f);
980 if (!faceRemoved(faceI) && f.size() < 3)
982 FatalErrorIn("polyTopoChange::compact(..)")
983 << "Created illegal face " << f
984 //<< " from face " << oldF
985 << " at position:" << faceI
986 << " when filtering removed points"
987 << abort(FatalError);
995 labelList localFaceMap(faces_.size(), -1);
998 forAll(faces_, faceI)
1000 if (!faceRemoved(faceI) && faceOwner_[faceI] >= 0)
1002 localFaceMap[faceI] = newFaceI++;
1005 nActiveFaces_ = newFaceI;
1007 forAll(faces_, faceI)
1009 if (!faceRemoved(faceI) && faceOwner_[faceI] < 0)
1012 localFaceMap[faceI] = newFaceI++;
1018 Pout<< "Faces : active:" << nActiveFaces_
1019 << " removed:" << faces_.size()-newFaceI << endl;
1023 reorderCompactFaces(newFaceI, localFaceMap);
1028 labelList localCellMap;
1033 // Construct cellCell addressing
1034 CompactListList<label> cellCells;
1035 makeCellCells(nActiveFaces_, cellCells);
1037 // Cell ordering (based on bandCompression). Handles removed cells.
1038 newCellI = getCellOrder(cellCells, localCellMap);
1042 // Compact out removed cells
1043 localCellMap.setSize(cellMap_.size());
1047 forAll(cellMap_, cellI)
1049 if (!cellRemoved(cellI))
1051 localCellMap[cellI] = newCellI++;
1058 Pout<< "Cells : active:" << newCellI
1059 << " removed:" << cellMap_.size()-newCellI << endl;
1062 // Renumber -if cells reordered or -if cells removed
1063 if (orderCells || (newCellI != cellMap_.size()))
1065 reorder(localCellMap, cellMap_);
1066 cellMap_.setCapacity(newCellI);
1067 renumberReverseMap(localCellMap, reverseCellMap_);
1069 reorder(localCellMap, cellZone_);
1070 cellZone_.setCapacity(newCellI);
1072 renumberKey(localCellMap, cellFromPoint_);
1073 renumberKey(localCellMap, cellFromEdge_);
1074 renumberKey(localCellMap, cellFromFace_);
1076 // Renumber owner/neighbour. Take into account if neighbour suddenly
1077 // gets lower cell than owner.
1078 forAll(faceOwner_, faceI)
1080 label own = faceOwner_[faceI];
1081 label nei = faceNeighbour_[faceI];
1086 faceOwner_[faceI] = localCellMap[own];
1090 // Update neighbour.
1091 faceNeighbour_[faceI] = localCellMap[nei];
1093 // Check if face needs reversing.
1096 faceNeighbour_[faceI] >= 0
1097 && faceNeighbour_[faceI] < faceOwner_[faceI]
1100 faces_[faceI] = faces_[faceI].reverseFace();
1101 Swap(faceOwner_[faceI], faceNeighbour_[faceI]);
1102 flipFaceFlux_[faceI] =
1104 flipFaceFlux_[faceI]
1108 faceZoneFlip_[faceI] =
1110 faceZoneFlip_[faceI]
1119 // Update neighbour.
1120 faceNeighbour_[faceI] = localCellMap[nei];
1126 // Reorder faces into upper-triangular and patch ordering
1128 // Create cells (packed storage)
1129 labelList cellFaces;
1130 labelList cellFaceOffsets;
1131 makeCells(nActiveFaces_, cellFaces, cellFaceOffsets);
1133 // Do upper triangular order and patch sorting
1134 labelList localFaceMap;
1147 reorderCompactFaces(localFaceMap.size(), localFaceMap);
1152 // Find faces to interpolate to create value for new face. Only used if
1153 // face was inflated from edge or point. Internal faces should only be
1154 // created from internal faces, external faces only from external faces
1155 // (and ideally the same patch)
1156 // Is bit problematic if there are no faces to select, i.e. in polyDualMesh
1157 // an internal face can be created from a boundary edge with no internal
1158 // faces connected to it.
1159 Foam::labelList Foam::polyTopoChange::selectFaces
1161 const primitiveMesh& mesh,
1162 const labelList& faceLabels,
1163 const bool internalFacesOnly
1168 forAll(faceLabels, i)
1170 label faceI = faceLabels[i];
1172 if (internalFacesOnly == mesh.isInternalFace(faceI))
1178 labelList collectedFaces;
1182 // Did not find any faces of the correct type so just use any old
1184 collectedFaces = faceLabels;
1188 collectedFaces.setSize(nFaces);
1192 forAll(faceLabels, i)
1194 label faceI = faceLabels[i];
1196 if (internalFacesOnly == mesh.isInternalFace(faceI))
1198 collectedFaces[nFaces++] = faceI;
1203 return collectedFaces;
1207 // Calculate pointMap per patch (so from patch point label to old patch point
1209 void Foam::polyTopoChange::calcPatchPointMap
1211 const List<Map<label> >& oldPatchMeshPointMaps,
1212 const polyBoundaryMesh& boundary,
1213 labelListList& patchPointMap
1216 patchPointMap.setSize(boundary.size());
1218 forAll(boundary, patchI)
1220 const labelList& meshPoints = boundary[patchI].meshPoints();
1222 const Map<label>& oldMeshPointMap = oldPatchMeshPointMaps[patchI];
1224 labelList& curPatchPointRnb = patchPointMap[patchI];
1226 curPatchPointRnb.setSize(meshPoints.size());
1228 forAll(meshPoints, i)
1230 if (meshPoints[i] < pointMap_.size())
1232 // Check if old point was part of same patch
1233 Map<label>::const_iterator ozmpmIter = oldMeshPointMap.find
1235 pointMap_[meshPoints[i]]
1238 if (ozmpmIter != oldMeshPointMap.end())
1240 curPatchPointRnb[i] = ozmpmIter();
1244 curPatchPointRnb[i] = -1;
1249 curPatchPointRnb[i] = -1;
1256 void Foam::polyTopoChange::calcFaceInflationMaps
1258 const polyMesh& mesh,
1259 List<objectMap>& facesFromPoints,
1260 List<objectMap>& facesFromEdges,
1261 List<objectMap>& facesFromFaces
1264 // Faces inflated from points
1265 // ~~~~~~~~~~~~~~~~~~~~~~~~~~
1267 facesFromPoints.setSize(faceFromPoint_.size());
1269 if (faceFromPoint_.size())
1271 label nFacesFromPoints = 0;
1273 // Collect all still existing faces connected to this point.
1274 forAllConstIter(Map<label>, faceFromPoint_, iter)
1276 label newFaceI = iter.key();
1278 if (region_[newFaceI] == -1)
1280 // Get internal faces using point on old mesh
1281 facesFromPoints[nFacesFromPoints++] = objectMap
1287 mesh.pointFaces()[iter()],
1294 // Get patch faces using point on old mesh
1295 facesFromPoints[nFacesFromPoints++] = objectMap
1301 mesh.pointFaces()[iter()],
1310 // Faces inflated from edges
1311 // ~~~~~~~~~~~~~~~~~~~~~~~~~
1313 facesFromEdges.setSize(faceFromEdge_.size());
1315 if (faceFromEdge_.size())
1317 label nFacesFromEdges = 0;
1319 // Collect all still existing faces connected to this edge.
1320 forAllConstIter(Map<label>, faceFromEdge_, iter)
1322 label newFaceI = iter.key();
1324 if (region_[newFaceI] == -1)
1326 // Get internal faces using edge on old mesh
1327 facesFromEdges[nFacesFromEdges++] = objectMap
1333 mesh.edgeFaces(iter()),
1340 // Get patch faces using edge on old mesh
1341 facesFromEdges[nFacesFromEdges++] = objectMap
1347 mesh.edgeFaces(iter()),
1356 // Faces from face merging
1357 // ~~~~~~~~~~~~~~~~~~~~~~~
1368 void Foam::polyTopoChange::calcCellInflationMaps
1370 const polyMesh& mesh,
1371 List<objectMap>& cellsFromPoints,
1372 List<objectMap>& cellsFromEdges,
1373 List<objectMap>& cellsFromFaces,
1374 List<objectMap>& cellsFromCells
1377 cellsFromPoints.setSize(cellFromPoint_.size());
1379 if (cellFromPoint_.size())
1381 label nCellsFromPoints = 0;
1383 // Collect all still existing faces connected to this point.
1384 forAllConstIter(Map<label>, cellFromPoint_, iter)
1386 cellsFromPoints[nCellsFromPoints++] = objectMap
1389 mesh.pointCells()[iter()]
1395 cellsFromEdges.setSize(cellFromEdge_.size());
1397 if (cellFromEdge_.size())
1399 label nCellsFromEdges = 0;
1401 // Collect all still existing faces connected to this point.
1402 forAllConstIter(Map<label>, cellFromEdge_, iter)
1404 cellsFromEdges[nCellsFromEdges++] = objectMap
1407 mesh.edgeCells()[iter()]
1413 cellsFromFaces.setSize(cellFromFace_.size());
1415 if (cellFromFace_.size())
1417 label nCellsFromFaces = 0;
1419 labelList cellsAroundFace(2, -1);
1421 // Collect all still existing faces connected to this point.
1422 forAllConstIter(Map<label>, cellFromFace_, iter)
1424 label oldFaceI = iter();
1426 if (mesh.isInternalFace(oldFaceI))
1428 cellsAroundFace[0] = mesh.faceOwner()[oldFaceI];
1429 cellsAroundFace[1] = mesh.faceNeighbour()[oldFaceI];
1433 cellsAroundFace[0] = mesh.faceOwner()[oldFaceI];
1434 cellsAroundFace[1] = -1;
1437 cellsFromFaces[nCellsFromFaces++] = objectMap
1446 // Cells from cell merging
1447 // ~~~~~~~~~~~~~~~~~~~~~~~
1458 void Foam::polyTopoChange::resetZones
1460 const polyMesh& mesh,
1462 labelListList& pointZoneMap,
1463 labelListList& faceZoneFaceMap,
1464 labelListList& cellZoneMap
1470 pointZoneMap.setSize(mesh.pointZones().size());
1472 const pointZoneMesh& pointZones = mesh.pointZones();
1474 // Count points per zone
1476 labelList nPoints(pointZones.size(), 0);
1478 forAllConstIter(Map<label>, pointZone_, iter)
1480 label zoneI = iter();
1482 if (zoneI < 0 || zoneI >= pointZones.size())
1486 "resetZones(const polyMesh&, polyMesh&, labelListList&"
1487 "labelListList&, labelListList&)"
1488 ) << "Illegal zoneID " << zoneI << " for point "
1489 << iter.key() << " coord " << mesh.points()[iter.key()]
1490 << abort(FatalError);
1495 // Distribute points per zone
1497 labelListList addressing(pointZones.size());
1498 forAll(addressing, zoneI)
1500 addressing[zoneI].setSize(nPoints[zoneI]);
1504 forAllConstIter(Map<label>, pointZone_, iter)
1506 label zoneI = iter();
1508 addressing[zoneI][nPoints[zoneI]++] = iter.key();
1510 // Sort the addressing
1511 forAll(addressing, zoneI)
1513 stableSort(addressing[zoneI]);
1516 // So now we both have old zones and the new addressing.
1517 // Invert the addressing to get pointZoneMap.
1518 forAll(addressing, zoneI)
1520 const pointZone& oldZone = pointZones[zoneI];
1521 const labelList& newZoneAddr = addressing[zoneI];
1523 labelList& curPzRnb = pointZoneMap[zoneI];
1524 curPzRnb.setSize(newZoneAddr.size());
1526 forAll(newZoneAddr, i)
1528 if (newZoneAddr[i] < pointMap_.size())
1530 curPzRnb[i] = oldZone.whichPoint(pointMap_[newZoneAddr[i]]);
1539 // Reset the addresing on the zone
1540 newMesh.pointZones().clearAddressing();
1541 forAll(newMesh.pointZones(), zoneI)
1545 Pout<< "pointZone:" << zoneI
1546 << " name:" << newMesh.pointZones()[zoneI].name()
1547 << " size:" << addressing[zoneI].size()
1551 newMesh.pointZones()[zoneI] = addressing[zoneI];
1559 faceZoneFaceMap.setSize(mesh.faceZones().size());
1561 const faceZoneMesh& faceZones = mesh.faceZones();
1563 labelList nFaces(faceZones.size(), 0);
1565 forAllConstIter(Map<label>, faceZone_, iter)
1567 label zoneI = iter();
1569 if (zoneI < 0 || zoneI >= faceZones.size())
1573 "resetZones(const polyMesh&, polyMesh&, labelListList&"
1574 "labelListList&, labelListList&)"
1575 ) << "Illegal zoneID " << zoneI << " for face "
1577 << abort(FatalError);
1582 labelListList addressing(faceZones.size());
1583 boolListList flipMode(faceZones.size());
1585 forAll(addressing, zoneI)
1587 addressing[zoneI].setSize(nFaces[zoneI]);
1588 flipMode[zoneI].setSize(nFaces[zoneI]);
1592 forAllConstIter(Map<label>, faceZone_, iter)
1594 label zoneI = iter();
1595 label faceI = iter.key();
1597 label index = nFaces[zoneI]++;
1599 addressing[zoneI][index] = faceI;
1600 flipMode[zoneI][index] = faceZoneFlip_[faceI];
1602 // Sort the addressing
1603 forAll(addressing, zoneI)
1606 sortedOrder(addressing[zoneI], newToOld);
1608 labelList newAddressing(addressing[zoneI].size());
1609 forAll(newAddressing, i)
1611 newAddressing[i] = addressing[zoneI][newToOld[i]];
1613 addressing[zoneI].transfer(newAddressing);
1616 boolList newFlipMode(flipMode[zoneI].size());
1617 forAll(newFlipMode, i)
1619 newFlipMode[i] = flipMode[zoneI][newToOld[i]];
1621 flipMode[zoneI].transfer(newFlipMode);
1625 // So now we both have old zones and the new addressing.
1626 // Invert the addressing to get faceZoneFaceMap.
1627 forAll(addressing, zoneI)
1629 const faceZone& oldZone = faceZones[zoneI];
1630 const labelList& newZoneAddr = addressing[zoneI];
1632 labelList& curFzFaceRnb = faceZoneFaceMap[zoneI];
1634 curFzFaceRnb.setSize(newZoneAddr.size());
1636 forAll(newZoneAddr, i)
1638 if (newZoneAddr[i] < faceMap_.size())
1641 oldZone.whichFace(faceMap_[newZoneAddr[i]]);
1645 curFzFaceRnb[i] = -1;
1651 // Reset the addresing on the zone
1652 newMesh.faceZones().clearAddressing();
1653 forAll(newMesh.faceZones(), zoneI)
1657 Pout<< "faceZone:" << zoneI
1658 << " name:" << newMesh.faceZones()[zoneI].name()
1659 << " size:" << addressing[zoneI].size()
1663 newMesh.faceZones()[zoneI].resetAddressing
1675 cellZoneMap.setSize(mesh.cellZones().size());
1677 const cellZoneMesh& cellZones = mesh.cellZones();
1679 labelList nCells(cellZones.size(), 0);
1681 forAll(cellZone_, cellI)
1683 label zoneI = cellZone_[cellI];
1685 if (zoneI >= cellZones.size())
1689 "resetZones(const polyMesh&, polyMesh&, labelListList&"
1690 "labelListList&, labelListList&)"
1691 ) << "Illegal zoneID " << zoneI << " for cell "
1692 << cellI << abort(FatalError);
1701 labelListList addressing(cellZones.size());
1702 forAll(addressing, zoneI)
1704 addressing[zoneI].setSize(nCells[zoneI]);
1708 forAll(cellZone_, cellI)
1710 label zoneI = cellZone_[cellI];
1714 addressing[zoneI][nCells[zoneI]++] = cellI;
1717 // Sort the addressing
1718 forAll(addressing, zoneI)
1720 stableSort(addressing[zoneI]);
1723 // So now we both have old zones and the new addressing.
1724 // Invert the addressing to get cellZoneMap.
1725 forAll(addressing, zoneI)
1727 const cellZone& oldZone = cellZones[zoneI];
1728 const labelList& newZoneAddr = addressing[zoneI];
1730 labelList& curCellRnb = cellZoneMap[zoneI];
1732 curCellRnb.setSize(newZoneAddr.size());
1734 forAll(newZoneAddr, i)
1736 if (newZoneAddr[i] < cellMap_.size())
1739 oldZone.whichCell(cellMap_[newZoneAddr[i]]);
1748 // Reset the addresing on the zone
1749 newMesh.cellZones().clearAddressing();
1750 forAll(newMesh.cellZones(), zoneI)
1754 Pout<< "cellZone:" << zoneI
1755 << " name:" << newMesh.cellZones()[zoneI].name()
1756 << " size:" << addressing[zoneI].size()
1760 newMesh.cellZones()[zoneI] = addressing[zoneI];
1766 void Foam::polyTopoChange::calcFaceZonePointMap
1768 const polyMesh& mesh,
1769 const List<Map<label> >& oldFaceZoneMeshPointMaps,
1770 labelListList& faceZonePointMap
1773 const faceZoneMesh& faceZones = mesh.faceZones();
1775 faceZonePointMap.setSize(faceZones.size());
1777 forAll(faceZones, zoneI)
1779 const faceZone& newZone = faceZones[zoneI];
1781 const labelList& newZoneMeshPoints = newZone().meshPoints();
1783 const Map<label>& oldZoneMeshPointMap = oldFaceZoneMeshPointMaps[zoneI];
1785 labelList& curFzPointRnb = faceZonePointMap[zoneI];
1787 curFzPointRnb.setSize(newZoneMeshPoints.size());
1789 forAll(newZoneMeshPoints, pointI)
1791 if (newZoneMeshPoints[pointI] < pointMap_.size())
1793 Map<label>::const_iterator ozmpmIter =
1794 oldZoneMeshPointMap.find
1796 pointMap_[newZoneMeshPoints[pointI]]
1799 if (ozmpmIter != oldZoneMeshPointMap.end())
1801 curFzPointRnb[pointI] = ozmpmIter();
1805 curFzPointRnb[pointI] = -1;
1810 curFzPointRnb[pointI] = -1;
1817 Foam::face Foam::polyTopoChange::rotateFace
1823 face newF(f.size());
1827 label fp1 = (fp + nPos) % f.size();
1841 void Foam::polyTopoChange::reorderCoupledFaces
1843 const bool syncParallel,
1844 const polyBoundaryMesh& boundary,
1845 const labelList& patchStarts,
1846 const labelList& patchSizes,
1847 const pointField& points
1850 // Mapping for faces (old to new). Extends over all mesh faces for
1851 // convenience (could be just the external faces)
1852 labelList oldToNew(identity(faces_.size()));
1854 // Rotation on new faces.
1855 labelList rotation(faces_.size(), 0);
1858 forAll(boundary, patchI)
1860 if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1862 boundary[patchI].initOrder
1878 // Receive and calculate ordering
1880 bool anyChanged = false;
1882 forAll(boundary, patchI)
1884 if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1886 labelList patchFaceMap(patchSizes[patchI], -1);
1887 labelList patchFaceRotation(patchSizes[patchI], 0);
1889 bool changed = boundary[patchI].order
1907 // Merge patch face reordering into mesh face reordering table
1908 label start = patchStarts[patchI];
1910 forAll(patchFaceMap, patchFaceI)
1912 oldToNew[patchFaceI + start] =
1913 start + patchFaceMap[patchFaceI];
1916 forAll(patchFaceRotation, patchFaceI)
1918 rotation[patchFaceI + start] =
1919 patchFaceRotation[patchFaceI];
1929 reduce(anyChanged, orOp<bool>());
1934 // Reorder faces according to oldToNew.
1935 reorderCompactFaces(oldToNew.size(), oldToNew);
1937 // Rotate faces (rotation is already in new face indices).
1938 forAll(rotation, faceI)
1940 if (rotation[faceI] != 0)
1942 faces_[faceI] = rotateFace(faces_[faceI], rotation[faceI]);
1949 void Foam::polyTopoChange::compactAndReorder
1951 const polyMesh& mesh,
1952 const bool syncParallel,
1953 const bool orderCells,
1954 const bool orderPoints,
1956 label& nInternalPoints,
1957 pointField& newPoints,
1958 labelList& patchSizes,
1959 labelList& patchStarts,
1960 List<objectMap>& pointsFromPoints,
1961 List<objectMap>& facesFromPoints,
1962 List<objectMap>& facesFromEdges,
1963 List<objectMap>& facesFromFaces,
1964 List<objectMap>& cellsFromPoints,
1965 List<objectMap>& cellsFromEdges,
1966 List<objectMap>& cellsFromFaces,
1967 List<objectMap>& cellsFromCells,
1968 List<Map<label> >& oldPatchMeshPointMaps,
1969 labelList& oldPatchNMeshPoints,
1970 labelList& oldPatchStarts,
1971 List<Map<label> >& oldFaceZoneMeshPointMaps
1974 if (mesh.boundaryMesh().size() != nPatches_)
1976 FatalErrorIn("polyTopoChange::compactAndReorder(..)")
1977 << "polyTopoChange was constructed with a mesh with "
1978 << nPatches_ << " patches." << endl
1979 << "The mesh now provided has a different number of patches "
1980 << mesh.boundaryMesh().size()
1981 << " which is illegal" << endl
1982 << abort(FatalError);
1985 // Remove any holes from points/faces/cells and sort faces.
1986 // Sets nActiveFaces_.
1987 compact(orderCells, orderPoints, nInternalPoints, patchSizes, patchStarts);
1989 // Transfer points to pointField. points_ are now cleared!
1990 // Only done since e.g. reorderCoupledFaces requires pointField.
1991 newPoints.transfer(points_);
1993 // Reorder any coupled faces
1997 mesh.boundaryMesh(),
2004 // Calculate inflation/merging maps
2005 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2006 // These are for the new face(/point/cell) the old faces whose value
2008 // averaged/summed to get the new value. These old faces come either from
2009 // merged old faces (face remove into other face),
2010 // the old edgeFaces (inflate from edge) or the old pointFaces (inflate
2011 // from point). As an additional complexity will use only internal faces
2012 // to create new value for internal face and vice versa only patch
2013 // faces to to create patch face value.
2015 // For point only point merging
2023 calcFaceInflationMaps
2031 calcCellInflationMaps
2040 // Clear inflation info
2042 faceFromPoint_.clearStorage();
2043 faceFromEdge_.clearStorage();
2045 cellFromPoint_.clearStorage();
2046 cellFromEdge_.clearStorage();
2047 cellFromFace_.clearStorage();
2051 const polyBoundaryMesh& boundary = mesh.boundaryMesh();
2053 // Grab patch mesh point maps
2054 oldPatchMeshPointMaps.setSize(boundary.size());
2055 oldPatchNMeshPoints.setSize(boundary.size());
2056 oldPatchStarts.setSize(boundary.size());
2058 forAll(boundary, patchI)
2060 // Copy old face zone mesh point maps
2061 oldPatchMeshPointMaps[patchI] = boundary[patchI].meshPointMap();
2062 oldPatchNMeshPoints[patchI] = boundary[patchI].meshPoints().size();
2063 oldPatchStarts[patchI] = boundary[patchI].start();
2066 // Grab old face zone mesh point maps.
2067 // These need to be saved before resetting the mesh and are used
2068 // later on to calculate the faceZone pointMaps.
2069 oldFaceZoneMeshPointMaps.setSize(mesh.faceZones().size());
2071 forAll(mesh.faceZones(), zoneI)
2073 const faceZone& oldZone = mesh.faceZones()[zoneI];
2075 oldFaceZoneMeshPointMaps[zoneI] = oldZone().meshPointMap();
2080 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
2082 // Construct from components
2083 Foam::polyTopoChange::polyTopoChange(const label nPatches, const bool strict)
2086 nPatches_(nPatches),
2089 reversePointMap_(0),
2113 // Construct from components
2114 Foam::polyTopoChange::polyTopoChange
2116 const polyMesh& mesh,
2124 reversePointMap_(0),
2149 identity(mesh.boundaryMesh().size()),
2150 identity(mesh.pointZones().size()),
2151 identity(mesh.faceZones().size()),
2152 identity(mesh.cellZones().size())
2157 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2159 void Foam::polyTopoChange::clear()
2161 points_.clearStorage();
2162 pointMap_.clearStorage();
2163 reversePointMap_.clearStorage();
2164 pointZone_.clearStorage();
2165 retiredPoints_.clearStorage();
2167 faces_.clearStorage();
2168 region_.clearStorage();
2169 faceOwner_.clearStorage();
2170 faceNeighbour_.clearStorage();
2171 faceMap_.clearStorage();
2172 reverseFaceMap_.clearStorage();
2173 faceFromPoint_.clearStorage();
2174 faceFromEdge_.clearStorage();
2175 flipFaceFlux_.clearStorage();
2176 faceZone_.clearStorage();
2177 faceZoneFlip_.clearStorage();
2180 cellMap_.clearStorage();
2181 reverseCellMap_.clearStorage();
2182 cellZone_.clearStorage();
2183 cellFromPoint_.clearStorage();
2184 cellFromEdge_.clearStorage();
2185 cellFromFace_.clearStorage();
2189 void Foam::polyTopoChange::addMesh
2191 const polyMesh& mesh,
2192 const labelList& patchMap,
2193 const labelList& pointZoneMap,
2194 const labelList& faceZoneMap,
2195 const labelList& cellZoneMap
2198 label maxRegion = nPatches_ - 1;
2201 maxRegion = max(maxRegion, patchMap[i]);
2203 nPatches_ = maxRegion + 1;
2208 const pointField& points = mesh.points();
2209 const pointZoneMesh& pointZones = mesh.pointZones();
2212 points_.setCapacity(points_.size() + points.size());
2213 pointMap_.setCapacity(pointMap_.size() + points.size());
2214 reversePointMap_.setCapacity(reversePointMap_.size() + points.size());
2215 pointZone_.resize(pointZone_.size() + points.size()/100);
2217 // Precalc offset zones
2218 labelList newZoneID(points.size(), -1);
2220 forAll(pointZones, zoneI)
2222 const labelList& pointLabels = pointZones[zoneI];
2224 forAll(pointLabels, j)
2226 newZoneID[pointLabels[j]] = pointZoneMap[zoneI];
2230 // Add points in mesh order
2231 for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
2245 const cellZoneMesh& cellZones = mesh.cellZones();
2249 // Note: polyMesh does not allow retired cells anymore. So allCells
2250 // always equals nCells
2251 label nAllCells = mesh.nCells();
2253 cellMap_.setCapacity(cellMap_.size() + nAllCells);
2254 reverseCellMap_.setCapacity(reverseCellMap_.size() + nAllCells);
2255 cellFromPoint_.resize(cellFromPoint_.size() + nAllCells/100);
2256 cellFromEdge_.resize(cellFromEdge_.size() + nAllCells/100);
2257 cellFromFace_.resize(cellFromFace_.size() + nAllCells/100);
2258 cellZone_.setCapacity(cellZone_.size() + nAllCells);
2261 // Precalc offset zones
2262 labelList newZoneID(nAllCells, -1);
2264 forAll(cellZones, zoneI)
2266 const labelList& cellLabels = cellZones[zoneI];
2268 forAll(cellLabels, j)
2270 label cellI = cellLabels[j];
2272 if (newZoneID[cellI] != -1)
2276 "polyTopoChange::addMesh"
2277 "(const polyMesh&, const labelList&,"
2278 "const labelList&, const labelList&,"
2280 ) << "Cell:" << cellI
2281 << " centre:" << mesh.cellCentres()[cellI]
2282 << " is in two zones:"
2283 << cellZones[newZoneID[cellI]].name()
2284 << " and " << cellZones[zoneI].name() << endl
2285 << " This is not supported."
2286 << " Continuing with first zone only." << endl;
2290 newZoneID[cellI] = cellZoneMap[zoneI];
2295 // Add cells in mesh order
2296 for (label cellI = 0; cellI < nAllCells; cellI++)
2298 // Add cell from cell
2299 addCell(-1, -1, -1, cellI, newZoneID[cellI]);
2305 const polyBoundaryMesh& patches = mesh.boundaryMesh();
2306 const faceList& faces = mesh.faces();
2307 const labelList& faceOwner = mesh.faceOwner();
2308 const labelList& faceNeighbour = mesh.faceNeighbour();
2309 const faceZoneMesh& faceZones = mesh.faceZones();
2312 label nAllFaces = mesh.faces().size();
2314 faces_.setCapacity(faces_.size() + nAllFaces);
2315 region_.setCapacity(region_.size() + nAllFaces);
2316 faceOwner_.setCapacity(faceOwner_.size() + nAllFaces);
2317 faceNeighbour_.setCapacity(faceNeighbour_.size() + nAllFaces);
2318 faceMap_.setCapacity(faceMap_.size() + nAllFaces);
2319 reverseFaceMap_.setCapacity(reverseFaceMap_.size() + nAllFaces);
2320 faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/100);
2321 faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/100);
2322 flipFaceFlux_.setCapacity(faces_.size() + nAllFaces);
2323 faceZone_.resize(faceZone_.size() + nAllFaces/100);
2324 faceZoneFlip_.setCapacity(faces_.size() + nAllFaces);
2327 // Precalc offset zones
2328 labelList newZoneID(nAllFaces, -1);
2329 boolList zoneFlip(nAllFaces, false);
2331 forAll(faceZones, zoneI)
2333 const labelList& faceLabels = faceZones[zoneI];
2334 const boolList& flipMap = faceZones[zoneI].flipMap();
2336 forAll(faceLabels, j)
2338 newZoneID[faceLabels[j]] = faceZoneMap[zoneI];
2339 zoneFlip[faceLabels[j]] = flipMap[j];
2343 // Add faces in mesh order
2345 // 1. Internal faces
2346 for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
2352 faceNeighbour[faceI],
2353 -1, // masterPointID
2355 faceI, // masterFaceID
2356 false, // flipFaceFlux
2358 newZoneID[faceI], // zoneID
2359 zoneFlip[faceI] // zoneFlip
2364 forAll(patches, patchI)
2366 const polyPatch& pp = patches[patchI];
2368 if (pp.start() != faces_.size())
2372 "polyTopoChange::polyTopoChange"
2373 "(const polyMesh& mesh, const bool strict)"
2375 << "Patch " << pp.name() << " starts at " << pp.start()
2377 << "Current face counter at " << faces_.size() << endl
2378 << "Are patches in incremental order?"
2379 << abort(FatalError);
2381 forAll(pp, patchFaceI)
2383 label faceI = pp.start() + patchFaceI;
2390 -1, // masterPointID
2392 faceI, // masterFaceID
2393 false, // flipFaceFlux
2394 patchMap[patchI], // patchID
2395 newZoneID[faceI], // zoneID
2396 zoneFlip[faceI] // zoneFlip
2404 void Foam::polyTopoChange::setCapacity
2406 const label nPoints,
2411 points_.setCapacity(nPoints);
2412 pointMap_.setCapacity(nPoints);
2413 reversePointMap_.setCapacity(nPoints);
2414 pointZone_.resize(pointZone_.size() + nPoints/100);
2416 faces_.setCapacity(nFaces);
2417 region_.setCapacity(nFaces);
2418 faceOwner_.setCapacity(nFaces);
2419 faceNeighbour_.setCapacity(nFaces);
2420 faceMap_.setCapacity(nFaces);
2421 reverseFaceMap_.setCapacity(nFaces);
2422 faceFromPoint_.resize(faceFromPoint_.size() + nFaces/100);
2423 faceFromEdge_.resize(faceFromEdge_.size() + nFaces/100);
2424 flipFaceFlux_.setCapacity(nFaces);
2425 faceZone_.resize(faceZone_.size() + nFaces/100);
2426 faceZoneFlip_.setCapacity(nFaces);
2428 cellMap_.setCapacity(nCells);
2429 reverseCellMap_.setCapacity(nCells);
2430 cellFromPoint_.resize(cellFromPoint_.size() + nCells/100);
2431 cellFromEdge_.resize(cellFromEdge_.size() + nCells/100);
2432 cellFromFace_.resize(cellFromFace_.size() + nCells/100);
2433 cellZone_.setCapacity(nCells);
2437 Foam::label Foam::polyTopoChange::setAction(const topoAction& action)
2439 if (isType<polyAddPoint>(action))
2441 const polyAddPoint& pap = refCast<const polyAddPoint>(action);
2446 pap.masterPointID(),
2451 else if (isType<polyModifyPoint>(action))
2453 const polyModifyPoint& pmp = refCast<const polyModifyPoint>(action);
2465 else if (isType<polyRemovePoint>(action))
2467 const polyRemovePoint& prp = refCast<const polyRemovePoint>(action);
2469 removePoint(prp.pointID(), prp.mergePointID());
2473 else if (isType<polyAddFace>(action))
2475 const polyAddFace& paf = refCast<const polyAddFace>(action);
2482 paf.masterPointID(),
2491 else if (isType<polyModifyFace>(action))
2493 const polyModifyFace& pmf = refCast<const polyModifyFace>(action);
2509 else if (isType<polyRemoveFace>(action))
2511 const polyRemoveFace& prf = refCast<const polyRemoveFace>(action);
2513 removeFace(prf.faceID(), prf.mergeFaceID());
2517 else if (isType<polyAddCell>(action))
2519 const polyAddCell& pac = refCast<const polyAddCell>(action);
2523 pac.masterPointID(),
2530 else if (isType<polyModifyCell>(action))
2532 const polyModifyCell& pmc = refCast<const polyModifyCell>(action);
2534 if (pmc.removeFromZone())
2536 modifyCell(pmc.cellID(), -1);
2540 modifyCell(pmc.cellID(), pmc.zoneID());
2545 else if (isType<polyRemoveCell>(action))
2547 const polyRemoveCell& prc = refCast<const polyRemoveCell>(action);
2549 removeCell(prc.cellID(), prc.mergeCellID());
2557 "label polyTopoChange::setAction(const topoAction& action)"
2558 ) << "Unknown type of topoChange: " << action.type()
2559 << abort(FatalError);
2561 // Dummy return to keep compiler happy
2567 Foam::label Foam::polyTopoChange::addPoint
2570 const label masterPointID,
2575 label pointI = points_.size();
2578 pointMap_.append(masterPointID);
2579 reversePointMap_.append(pointI);
2583 pointZone_.insert(pointI, zoneID);
2588 retiredPoints_.insert(pointI);
2595 void Foam::polyTopoChange::modifyPoint
2599 const label newZoneID,
2603 if (pointI < 0 || pointI >= points_.size())
2607 "polyTopoChange::modifyPoint(const label, const point&)"
2608 ) << "illegal point label " << pointI << endl
2609 << "Valid point labels are 0 .. " << points_.size()-1
2610 << abort(FatalError);
2612 if (pointRemoved(pointI) || pointMap_[pointI] == -1)
2616 "polyTopoChange::modifyPoint(const label, const point&)"
2617 ) << "point " << pointI << " already marked for removal"
2618 << abort(FatalError);
2620 points_[pointI] = pt;
2622 Map<label>::iterator pointFnd = pointZone_.find(pointI);
2624 if (pointFnd != pointZone_.end())
2628 pointFnd() = newZoneID;
2632 pointZone_.erase(pointFnd);
2635 else if (newZoneID >= 0)
2637 pointZone_.insert(pointI, newZoneID);
2642 retiredPoints_.erase(pointI);
2646 retiredPoints_.insert(pointI);
2651 void Foam::polyTopoChange::movePoints(const pointField& newPoints)
2653 if (newPoints.size() != points_.size())
2655 FatalErrorIn("polyTopoChange::movePoints(const pointField&)")
2656 << "illegal pointField size." << endl
2657 << "Size:" << newPoints.size() << endl
2658 << "Points in mesh:" << points_.size()
2659 << abort(FatalError);
2662 forAll(points_, pointI)
2664 points_[pointI] = newPoints[pointI];
2669 void Foam::polyTopoChange::removePoint
2672 const label mergePointI
2675 if (pointI < 0 || pointI >= points_.size())
2677 FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2678 << "illegal point label " << pointI << endl
2679 << "Valid point labels are 0 .. " << points_.size()-1
2680 << abort(FatalError);
2686 && (pointRemoved(pointI) || pointMap_[pointI] == -1)
2689 FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2690 << "point " << pointI << " already marked for removal" << nl
2691 << "Point:" << points_[pointI] << " pointMap:" << pointMap_[pointI]
2692 << abort(FatalError);
2695 if (pointI == mergePointI)
2697 FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2698 << "Cannot remove/merge point " << pointI << " onto itself."
2699 << abort(FatalError);
2702 points_[pointI] = greatPoint;
2703 pointMap_[pointI] = -1;
2704 if (mergePointI >= 0)
2706 reversePointMap_[pointI] = -mergePointI-2;
2710 reversePointMap_[pointI] = -1;
2712 pointZone_.erase(pointI);
2713 retiredPoints_.erase(pointI);
2717 Foam::label Foam::polyTopoChange::addFace
2722 const label masterPointID,
2723 const label masterEdgeID,
2724 const label masterFaceID,
2725 const bool flipFaceFlux,
2726 const label patchID,
2734 checkFace(f, -1, own, nei, patchID, zoneID);
2737 label faceI = faces_.size();
2740 region_.append(patchID);
2741 faceOwner_.append(own);
2742 faceNeighbour_.append(nei);
2744 if (masterPointID >= 0)
2746 faceMap_.append(-1);
2747 faceFromPoint_.insert(faceI, masterPointID);
2749 else if (masterEdgeID >= 0)
2751 faceMap_.append(-1);
2752 faceFromEdge_.insert(faceI, masterEdgeID);
2754 else if (masterFaceID >= 0)
2756 faceMap_.append(masterFaceID);
2760 // Allow inflate-from-nothing?
2761 //FatalErrorIn("polyTopoChange::addFace")
2762 // << "Need to specify a master point, edge or face"
2763 // << "face:" << f << " own:" << own << " nei:" << nei
2764 // << abort(FatalError);
2765 faceMap_.append(-1);
2767 reverseFaceMap_.append(faceI);
2769 flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2773 faceZone_.insert(faceI, zoneID);
2775 faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2781 void Foam::polyTopoChange::modifyFace
2787 const bool flipFaceFlux,
2788 const label patchID,
2796 checkFace(f, faceI, own, nei, patchID, zoneID);
2800 faceOwner_[faceI] = own;
2801 faceNeighbour_[faceI] = nei;
2802 region_[faceI] = patchID;
2804 flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2806 Map<label>::iterator faceFnd = faceZone_.find(faceI);
2808 if (faceFnd != faceZone_.end())
2816 faceZone_.erase(faceFnd);
2819 else if (zoneID >= 0)
2821 faceZone_.insert(faceI, zoneID);
2823 faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2827 void Foam::polyTopoChange::removeFace(const label faceI, const label mergeFaceI)
2829 if (faceI < 0 || faceI >= faces_.size())
2831 FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2832 << "illegal face label " << faceI << endl
2833 << "Valid face labels are 0 .. " << faces_.size()-1
2834 << abort(FatalError);
2840 && (faceRemoved(faceI) || faceMap_[faceI] == -1)
2843 FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2845 << " already marked for removal"
2846 << abort(FatalError);
2849 faces_[faceI].setSize(0);
2850 region_[faceI] = -1;
2851 faceOwner_[faceI] = -1;
2852 faceNeighbour_[faceI] = -1;
2853 faceMap_[faceI] = -1;
2854 if (mergeFaceI >= 0)
2856 reverseFaceMap_[faceI] = -mergeFaceI-2;
2860 reverseFaceMap_[faceI] = -1;
2862 faceFromEdge_.erase(faceI);
2863 faceFromPoint_.erase(faceI);
2864 flipFaceFlux_[faceI] = 0;
2865 faceZone_.erase(faceI);
2866 faceZoneFlip_[faceI] = 0;
2870 Foam::label Foam::polyTopoChange::addCell
2872 const label masterPointID,
2873 const label masterEdgeID,
2874 const label masterFaceID,
2875 const label masterCellID,
2879 label cellI = cellMap_.size();
2881 if (masterPointID >= 0)
2883 cellMap_.append(-1);
2884 cellFromPoint_.insert(cellI, masterPointID);
2886 else if (masterEdgeID >= 0)
2888 cellMap_.append(-1);
2889 cellFromEdge_.insert(cellI, masterEdgeID);
2891 else if (masterFaceID >= 0)
2893 cellMap_.append(-1);
2894 cellFromFace_.insert(cellI, masterFaceID);
2898 cellMap_.append(masterCellID);
2900 reverseCellMap_.append(cellI);
2901 cellZone_.append(zoneID);
2907 void Foam::polyTopoChange::modifyCell
2913 cellZone_[cellI] = zoneID;
2917 void Foam::polyTopoChange::removeCell(const label cellI, const label mergeCellI)
2919 if (cellI < 0 || cellI >= cellMap_.size())
2921 FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2922 << "illegal cell label " << cellI << endl
2923 << "Valid cell labels are 0 .. " << cellMap_.size()-1
2924 << abort(FatalError);
2927 if (strict_ && cellMap_[cellI] == -2)
2929 FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2931 << " already marked for removal"
2932 << abort(FatalError);
2935 cellMap_[cellI] = -2;
2936 if (mergeCellI >= 0)
2938 reverseCellMap_[cellI] = -mergeCellI-2;
2942 reverseCellMap_[cellI] = -1;
2944 cellFromPoint_.erase(cellI);
2945 cellFromEdge_.erase(cellI);
2946 cellFromFace_.erase(cellI);
2947 cellZone_[cellI] = -1;
2951 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
2955 const bool syncParallel,
2956 const bool orderCells,
2957 const bool orderPoints
2962 Pout<< "polyTopoChange::changeMesh"
2963 << "(polyMesh&, const bool, const bool, const bool, const bool)"
2969 Pout<< "Old mesh:" << nl;
2970 writeMeshStats(mesh, Pout);
2974 pointField newPoints;
2975 // number of internal points
2976 label nInternalPoints;
2978 labelList patchSizes;
2979 labelList patchStarts;
2981 List<objectMap> pointsFromPoints;
2982 List<objectMap> facesFromPoints;
2983 List<objectMap> facesFromEdges;
2984 List<objectMap> facesFromFaces;
2985 List<objectMap> cellsFromPoints;
2986 List<objectMap> cellsFromEdges;
2987 List<objectMap> cellsFromFaces;
2988 List<objectMap> cellsFromCells;
2990 List<Map<label> > oldPatchMeshPointMaps;
2991 labelList oldPatchNMeshPoints;
2992 labelList oldPatchStarts;
2993 List<Map<label> > oldFaceZoneMeshPointMaps;
2995 // Compact, reorder patch faces and calculate mesh/patch maps.
3015 oldPatchMeshPointMaps,
3016 oldPatchNMeshPoints,
3018 oldFaceZoneMeshPointMaps
3021 const label nOldPoints(mesh.nPoints());
3022 const label nOldFaces(mesh.nFaces());
3023 const label nOldCells(mesh.nCells());
3028 // This will invalidate any addressing so better make sure you have
3029 // all the information you need!!!
3033 // Keep (renumbered) mesh points, store new points in map for inflation
3034 // (appended points (i.e. from nowhere) get value zero)
3035 pointField renumberedMeshPoints(newPoints.size());
3037 forAll(pointMap_, newPointI)
3039 label oldPointI = pointMap_[newPointI];
3043 renumberedMeshPoints[newPointI] = mesh.points()[oldPointI];
3047 renumberedMeshPoints[newPointI] = vector::zero;
3051 mesh.resetPrimitives
3053 xferMove(renumberedMeshPoints),
3056 faceNeighbour_.xfer(),
3062 mesh.changing(true);
3067 mesh.resetPrimitives
3069 xferMove(newPoints),
3072 faceNeighbour_.xfer(),
3077 mesh.changing(true);
3080 // Clear out primitives
3082 retiredPoints_.clearStorage();
3083 region_.clearStorage();
3089 // Some stats on changes
3090 label nAdd, nInflate, nMerge, nRemove;
3091 countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3093 << " added(from point):" << nAdd
3094 << " added(from nothing):" << nInflate
3095 << " merged(into other point):" << nMerge
3096 << " removed:" << nRemove
3099 countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3101 << " added(from face):" << nAdd
3102 << " added(inflated):" << nInflate
3103 << " merged(into other face):" << nMerge
3104 << " removed:" << nRemove
3107 countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3109 << " added(from cell):" << nAdd
3110 << " added(inflated):" << nInflate
3111 << " merged(into other cell):" << nMerge
3112 << " removed:" << nRemove
3119 Pout<< "New mesh:" << nl;
3120 writeMeshStats(mesh, Pout);
3127 // Inverse of point/face/cell zone addressing.
3128 // For every preserved point/face/cells in zone give the old position.
3129 // For added points, the index is set to -1
3130 labelListList pointZoneMap(mesh.pointZones().size());
3131 labelListList faceZoneFaceMap(mesh.faceZones().size());
3132 labelListList cellZoneMap(mesh.cellZones().size());
3134 resetZones(mesh, mesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3138 pointZone_.clearStorage();
3139 faceZone_.clearStorage();
3140 faceZoneFlip_.clearStorage();
3141 cellZone_.clearStorage();
3145 // Patch point renumbering
3146 // For every preserved point on a patch give the old position.
3147 // For added points, the index is set to -1
3148 labelListList patchPointMap(mesh.boundaryMesh().size());
3151 oldPatchMeshPointMaps,
3152 mesh.boundaryMesh(),
3156 // Create the face zone mesh point renumbering
3157 labelListList faceZonePointMap(mesh.faceZones().size());
3158 calcFaceZonePointMap(mesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3160 labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3162 return autoPtr<mapPolyMesh>
3199 newPoints, // if empty signals no inflation.
3201 oldPatchNMeshPoints,
3202 true // steal storage.
3206 // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3211 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
3213 autoPtr<fvMesh>& newMeshPtr,
3216 const bool syncParallel,
3217 const bool orderCells,
3218 const bool orderPoints
3223 Pout<< "polyTopoChange::changeMesh"
3224 << "(autoPtr<fvMesh>&, const IOobject&, const fvMesh&"
3225 << ", const bool, const bool, const bool)"
3231 Pout<< "Old mesh:" << nl;
3232 writeMeshStats(mesh, Pout);
3236 pointField newPoints;
3237 // number of internal points
3238 label nInternalPoints;
3240 labelList patchSizes;
3241 labelList patchStarts;
3243 List<objectMap> pointsFromPoints;
3244 List<objectMap> facesFromPoints;
3245 List<objectMap> facesFromEdges;
3246 List<objectMap> facesFromFaces;
3247 List<objectMap> cellsFromPoints;
3248 List<objectMap> cellsFromEdges;
3249 List<objectMap> cellsFromFaces;
3250 List<objectMap> cellsFromCells;
3253 List<Map<label> > oldPatchMeshPointMaps;
3254 labelList oldPatchNMeshPoints;
3255 labelList oldPatchStarts;
3256 List<Map<label> > oldFaceZoneMeshPointMaps;
3258 // Compact, reorder patch faces and calculate mesh/patch maps.
3278 oldPatchMeshPointMaps,
3279 oldPatchNMeshPoints,
3281 oldFaceZoneMeshPointMaps
3284 const label nOldPoints(mesh.nPoints());
3285 const label nOldFaces(mesh.nFaces());
3286 const label nOldCells(mesh.nCells());
3297 xferMove(newPoints),
3300 faceNeighbour_.xfer()
3303 fvMesh& newMesh = newMeshPtr();
3305 // Clear out primitives
3307 retiredPoints_.clearStorage();
3308 region_.clearStorage();
3314 // Some stats on changes
3315 label nAdd, nInflate, nMerge, nRemove;
3316 countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3318 << " added(from point):" << nAdd
3319 << " added(from nothing):" << nInflate
3320 << " merged(into other point):" << nMerge
3321 << " removed:" << nRemove
3324 countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3326 << " added(from face):" << nAdd
3327 << " added(inflated):" << nInflate
3328 << " merged(into other face):" << nMerge
3329 << " removed:" << nRemove
3332 countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3334 << " added(from cell):" << nAdd
3335 << " added(inflated):" << nInflate
3336 << " merged(into other cell):" << nMerge
3337 << " removed:" << nRemove
3344 const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
3346 List<polyPatch*> newBoundary(oldPatches.size());
3348 forAll(oldPatches, patchI)
3350 newBoundary[patchI] = oldPatches[patchI].clone
3352 newMesh.boundaryMesh(),
3358 newMesh.addFvPatches(newBoundary);
3365 // Start off from empty zones.
3366 const pointZoneMesh& oldPointZones = mesh.pointZones();
3367 List<pointZone*> pZonePtrs(oldPointZones.size());
3369 forAll(oldPointZones, i)
3371 pZonePtrs[i] = new pointZone
3373 oldPointZones[i].name(),
3376 newMesh.pointZones()
3381 const faceZoneMesh& oldFaceZones = mesh.faceZones();
3382 List<faceZone*> fZonePtrs(oldFaceZones.size());
3384 forAll(oldFaceZones, i)
3386 fZonePtrs[i] = new faceZone
3388 oldFaceZones[i].name(),
3397 const cellZoneMesh& oldCellZones = mesh.cellZones();
3398 List<cellZone*> cZonePtrs(oldCellZones.size());
3400 forAll(oldCellZones, i)
3402 cZonePtrs[i] = new cellZone
3404 oldCellZones[i].name(),
3412 newMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
3414 // Inverse of point/face/cell zone addressing.
3415 // For every preserved point/face/cells in zone give the old position.
3416 // For added points, the index is set to -1
3417 labelListList pointZoneMap(mesh.pointZones().size());
3418 labelListList faceZoneFaceMap(mesh.faceZones().size());
3419 labelListList cellZoneMap(mesh.cellZones().size());
3421 resetZones(mesh, newMesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3425 pointZone_.clearStorage();
3426 faceZone_.clearStorage();
3427 faceZoneFlip_.clearStorage();
3428 cellZone_.clearStorage();
3431 // Patch point renumbering
3432 // For every preserved point on a patch give the old position.
3433 // For added points, the index is set to -1
3434 labelListList patchPointMap(newMesh.boundaryMesh().size());
3437 oldPatchMeshPointMaps,
3438 newMesh.boundaryMesh(),
3442 // Create the face zone mesh point renumbering
3443 labelListList faceZonePointMap(newMesh.faceZones().size());
3444 calcFaceZonePointMap(newMesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3448 Pout<< "New mesh:" << nl;
3449 writeMeshStats(mesh, Pout);
3452 labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3454 return autoPtr<mapPolyMesh>
3491 newPoints, // if empty signals no inflation.
3493 oldPatchNMeshPoints,
3494 true // steal storage.
3498 // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3503 // ************************************************************************* //