00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "polyTopoChange.H"
00027 #include <OpenFOAM/SortableList.H>
00028 #include <OpenFOAM/polyMesh.H>
00029 #include <dynamicMesh/polyAddPoint.H>
00030 #include <dynamicMesh/polyModifyPoint.H>
00031 #include <dynamicMesh/polyRemovePoint.H>
00032 #include <dynamicMesh/polyAddFace.H>
00033 #include <dynamicMesh/polyModifyFace.H>
00034 #include <dynamicMesh/polyRemoveFace.H>
00035 #include <dynamicMesh/polyAddCell.H>
00036 #include <dynamicMesh/polyModifyCell.H>
00037 #include <dynamicMesh/polyRemoveCell.H>
00038 #include <OpenFOAM/objectMap.H>
00039 #include <OpenFOAM/processorPolyPatch.H>
00040 #include <finiteVolume/fvMesh.H>
00041
00042
00043
00044 namespace Foam
00045 {
00046 defineTypeNameAndDebug(polyTopoChange, 0);
00047 }
00048
00049
00050 const Foam::point Foam::polyTopoChange::greatPoint
00051 (
00052 GREAT,
00053 GREAT,
00054 GREAT
00055 );
00056
00057
00058
00059
00060
00061 void Foam::polyTopoChange::renumberReverseMap
00062 (
00063 const labelList& map,
00064 DynamicList<label>& elems
00065 )
00066 {
00067 forAll(elems, elemI)
00068 {
00069 label val = elems[elemI];
00070
00071 if (val >= 0)
00072 {
00073 elems[elemI] = map[val];
00074 }
00075 else if (val < -1)
00076 {
00077 label mergedVal = -val-2;
00078 elems[elemI] = -map[mergedVal]-2;
00079 }
00080 }
00081 }
00082
00083
00084 void Foam::polyTopoChange::renumber
00085 (
00086 const labelList& map,
00087 labelHashSet& elems
00088 )
00089 {
00090 labelHashSet newElems(elems.size());
00091
00092 forAllConstIter(labelHashSet, elems, iter)
00093 {
00094 label newElem = map[iter.key()];
00095
00096 if (newElem >= 0)
00097 {
00098 newElems.insert(newElem);
00099 }
00100 }
00101
00102 elems.transfer(newElems);
00103 }
00104
00105
00106
00107 void Foam::polyTopoChange::renumberCompact
00108 (
00109 const labelList& map,
00110 labelList& elems
00111 )
00112 {
00113 label newElemI = 0;
00114
00115 forAll(elems, elemI)
00116 {
00117 label newVal = map[elems[elemI]];
00118
00119 if (newVal != -1)
00120 {
00121 elems[newElemI++] = newVal;
00122 }
00123 }
00124 elems.setSize(newElemI);
00125 }
00126
00127
00128 void Foam::polyTopoChange::countMap
00129 (
00130 const labelList& map,
00131 const labelList& reverseMap,
00132 label& nAdd,
00133 label& nInflate,
00134 label& nMerge,
00135 label& nRemove
00136 )
00137 {
00138 nAdd = 0;
00139 nInflate = 0;
00140 nMerge = 0;
00141 nRemove = 0;
00142
00143 forAll(map, newCellI)
00144 {
00145 label oldCellI = map[newCellI];
00146
00147 if (oldCellI >= 0)
00148 {
00149 if (reverseMap[oldCellI] == newCellI)
00150 {
00151
00152 }
00153 else
00154 {
00155
00156 nAdd++;
00157 }
00158 }
00159 else if (oldCellI == -1)
00160 {
00161
00162 nInflate++;
00163 }
00164 else
00165 {
00166 FatalErrorIn("countMap") << "old:" << oldCellI
00167 << " new:" << newCellI << abort(FatalError);
00168 }
00169 }
00170
00171 forAll(reverseMap, oldCellI)
00172 {
00173 label newCellI = reverseMap[oldCellI];
00174
00175 if (newCellI >= 0)
00176 {
00177
00178 }
00179 else if (newCellI == -1)
00180 {
00181
00182 nRemove++;
00183 }
00184 else
00185 {
00186
00187 nMerge++;
00188 }
00189 }
00190 }
00191
00192
00193 Foam::labelHashSet Foam::polyTopoChange::getSetIndices
00194 (
00195 const PackedBoolList& lst
00196 )
00197 {
00198 labelHashSet values(lst.count());
00199 forAll(lst, i)
00200 {
00201 if (lst[i])
00202 {
00203 values.insert(i);
00204 }
00205 }
00206 return values;
00207 }
00208
00209
00210 void Foam::polyTopoChange::writeMeshStats(const polyMesh& mesh, Ostream& os)
00211 {
00212 const polyBoundaryMesh& patches = mesh.boundaryMesh();
00213
00214 labelList patchSizes(patches.size());
00215 labelList patchStarts(patches.size());
00216 forAll(patches, patchI)
00217 {
00218 patchSizes[patchI] = patches[patchI].size();
00219 patchStarts[patchI] = patches[patchI].start();
00220 }
00221
00222 os << " Points : " << mesh.nPoints() << nl
00223 << " Faces : " << mesh.nFaces() << nl
00224 << " Cells : " << mesh.nCells() << nl
00225 << " PatchSizes : " << patchSizes << nl
00226 << " PatchStarts : " << patchStarts << nl
00227 << endl;
00228 }
00229
00230
00231 void Foam::polyTopoChange::getMergeSets
00232 (
00233 const labelList& reverseCellMap,
00234 const labelList& cellMap,
00235 List<objectMap>& cellsFromCells
00236 )
00237 {
00238
00239 labelList nMerged(cellMap.size(), 1);
00240
00241 forAll(reverseCellMap, oldCellI)
00242 {
00243 label newCellI = reverseCellMap[oldCellI];
00244
00245 if (newCellI < -1)
00246 {
00247 label mergeCellI = -newCellI-2;
00248
00249 nMerged[mergeCellI]++;
00250 }
00251 }
00252
00253
00254 labelList cellToMergeSet(cellMap.size(), -1);
00255
00256 label nSets = 0;
00257
00258 forAll(nMerged, cellI)
00259 {
00260 if (nMerged[cellI] > 1)
00261 {
00262 cellToMergeSet[cellI] = nSets++;
00263 }
00264 }
00265
00266
00267
00268
00269
00270
00271
00272 cellsFromCells.setSize(nSets);
00273
00274 forAll(reverseCellMap, oldCellI)
00275 {
00276 label newCellI = reverseCellMap[oldCellI];
00277
00278 if (newCellI < -1)
00279 {
00280 label mergeCellI = -newCellI-2;
00281
00282
00283
00284 label setI = cellToMergeSet[mergeCellI];
00285
00286 objectMap& mergeSet = cellsFromCells[setI];
00287
00288 if (mergeSet.masterObjects().empty())
00289 {
00290
00291
00292 mergeSet.index() = mergeCellI;
00293 mergeSet.masterObjects().setSize(nMerged[mergeCellI]);
00294
00295
00296 mergeSet.masterObjects()[0] = cellMap[mergeCellI];
00297
00298
00299 mergeSet.masterObjects()[1] = oldCellI;
00300
00301 nMerged[mergeCellI] = 2;
00302 }
00303 else
00304 {
00305 mergeSet.masterObjects()[nMerged[mergeCellI]++] = oldCellI;
00306 }
00307 }
00308 }
00309 }
00310
00311
00312 void Foam::polyTopoChange::checkFace
00313 (
00314 const face& f,
00315 const label faceI,
00316 const label own,
00317 const label nei,
00318 const label patchI,
00319 const label zoneI
00320 ) const
00321 {
00322 if (nei == -1)
00323 {
00324 if (own == -1 && zoneI != -1)
00325 {
00326
00327 }
00328 else if (patchI == -1 || patchI >= nPatches_)
00329 {
00330 FatalErrorIn
00331 (
00332 "polyTopoChange::checkFace(const face&, const label"
00333 ", const label, const label, const label)"
00334 ) << "Face has no neighbour (so external) but does not have"
00335 << " a valid patch" << nl
00336 << "f:" << f
00337 << " faceI(-1 if added face):" << faceI
00338 << " own:" << own << " nei:" << nei
00339 << " patchI:" << patchI << abort(FatalError);
00340 }
00341 }
00342 else
00343 {
00344 if (patchI != -1)
00345 {
00346 FatalErrorIn
00347 (
00348 "polyTopoChange::checkFace(const face&, const label"
00349 ", const label, const label, const label)"
00350 ) << "Cannot both have valid patchI and neighbour" << nl
00351 << "f:" << f
00352 << " faceI(-1 if added face):" << faceI
00353 << " own:" << own << " nei:" << nei
00354 << " patchI:" << patchI << abort(FatalError);
00355 }
00356
00357 if (nei <= own)
00358 {
00359 FatalErrorIn
00360 (
00361 "polyTopoChange::checkFace(const face&, const label"
00362 ", const label, const label, const label)"
00363 ) << "Owner cell label should be less than neighbour cell label"
00364 << nl
00365 << "f:" << f
00366 << " faceI(-1 if added face):" << faceI
00367 << " own:" << own << " nei:" << nei
00368 << " patchI:" << patchI << abort(FatalError);
00369 }
00370 }
00371
00372 if (f.size() < 3 || findIndex(f, -1) != -1)
00373 {
00374 FatalErrorIn
00375 (
00376 "polyTopoChange::checkFace(const face&, const label"
00377 ", const label, const label, const label)"
00378 ) << "Illegal vertices in face"
00379 << nl
00380 << "f:" << f
00381 << " faceI(-1 if added face):" << faceI
00382 << " own:" << own << " nei:" << nei
00383 << " patchI:" << patchI << abort(FatalError);
00384 }
00385 if (faceI >= 0 && faceI < faces_.size() && faceRemoved(faceI))
00386 {
00387 FatalErrorIn
00388 (
00389 "polyTopoChange::checkFace(const face&, const label"
00390 ", const label, const label, const label)"
00391 ) << "Face already marked for removal"
00392 << nl
00393 << "f:" << f
00394 << " faceI(-1 if added face):" << faceI
00395 << " own:" << own << " nei:" << nei
00396 << " patchI:" << patchI << abort(FatalError);
00397 }
00398 forAll(f, fp)
00399 {
00400 if (f[fp] < points_.size() && pointRemoved(f[fp]))
00401 {
00402 FatalErrorIn
00403 (
00404 "polyTopoChange::checkFace(const face&, const label"
00405 ", const label, const label, const label)"
00406 ) << "Face uses removed vertices"
00407 << nl
00408 << "f:" << f
00409 << " faceI(-1 if added face):" << faceI
00410 << " own:" << own << " nei:" << nei
00411 << " patchI:" << patchI << abort(FatalError);
00412 }
00413 }
00414 }
00415
00416
00417 void Foam::polyTopoChange::makeCells
00418 (
00419 const label nActiveFaces,
00420 labelList& cellFaces,
00421 labelList& cellFaceOffsets
00422 ) const
00423 {
00424 cellFaces.setSize(2*nActiveFaces);
00425 cellFaceOffsets.setSize(cellMap_.size() + 1);
00426
00427
00428 labelList nNbrs(cellMap_.size(), 0);
00429
00430
00431
00432 for (label faceI = 0; faceI < nActiveFaces; faceI++)
00433 {
00434 nNbrs[faceOwner_[faceI]]++;
00435 }
00436 for (label faceI = 0; faceI < nActiveFaces; faceI++)
00437 {
00438 if (faceNeighbour_[faceI] >= 0)
00439 {
00440 nNbrs[faceNeighbour_[faceI]]++;
00441 }
00442 }
00443
00444
00445
00446 cellFaceOffsets[0] = 0;
00447 forAll(nNbrs, cellI)
00448 {
00449 cellFaceOffsets[cellI+1] = cellFaceOffsets[cellI] + nNbrs[cellI];
00450 }
00451
00452
00453
00454
00455 nNbrs = 0;
00456
00457 for (label faceI = 0; faceI < nActiveFaces; faceI++)
00458 {
00459 label cellI = faceOwner_[faceI];
00460
00461 cellFaces[cellFaceOffsets[cellI] + nNbrs[cellI]++] = faceI;
00462 }
00463
00464 for (label faceI = 0; faceI < nActiveFaces; faceI++)
00465 {
00466 label cellI = faceNeighbour_[faceI];
00467
00468 if (cellI >= 0)
00469 {
00470 cellFaces[cellFaceOffsets[cellI] + nNbrs[cellI]++] = faceI;
00471 }
00472 }
00473
00474
00475 cellFaces.setSize(cellFaceOffsets[cellMap_.size()]);
00476 }
00477
00478
00479
00480
00481 void Foam::polyTopoChange::makeCellCells
00482 (
00483 const label nActiveFaces,
00484 CompactListList<label>& cellCells
00485 ) const
00486 {
00487
00488 labelList nNbrs(cellMap_.size(), 0);
00489
00490
00491 label nCellCells = 0;
00492
00493
00494
00495 for (label faceI = 0; faceI < nActiveFaces; faceI++)
00496 {
00497 if (faceNeighbour_[faceI] >= 0)
00498 {
00499 nNbrs[faceOwner_[faceI]]++;
00500 nNbrs[faceNeighbour_[faceI]]++;
00501 nCellCells += 2;
00502 }
00503 }
00504
00505 cellCells.setSize(cellMap_.size(), nCellCells);
00506
00507
00508
00509 labelList& offsets = cellCells.offsets();
00510
00511 label sumSize = 0;
00512 forAll(nNbrs, cellI)
00513 {
00514 sumSize += nNbrs[cellI];
00515 offsets[cellI] = sumSize;
00516 }
00517
00518
00519
00520
00521 nNbrs = 0;
00522
00523 for (label faceI = 0; faceI < nActiveFaces; faceI++)
00524 {
00525 label nei = faceNeighbour_[faceI];
00526
00527 if (nei >= 0)
00528 {
00529 label own = faceOwner_[faceI];
00530 cellCells.m()[cellCells.index(own, nNbrs[own]++)] = nei;
00531 cellCells.m()[cellCells.index(nei, nNbrs[nei]++)] = own;
00532 }
00533 }
00534 }
00535
00536
00537
00538
00539 Foam::label Foam::polyTopoChange::getCellOrder
00540 (
00541 const CompactListList<label>& cellCellAddressing,
00542 labelList& oldToNew
00543 ) const
00544 {
00545 const labelList& offsets = cellCellAddressing.offsets();
00546
00547 labelList newOrder(cellCellAddressing.size());
00548
00549
00550 SLList<label> nextCell;
00551
00552
00553 PackedBoolList visited(cellCellAddressing.size());
00554
00555 label cellInOrder = 0;
00556
00557
00558
00559 forAll (visited, cellI)
00560 {
00561
00562 if (!cellRemoved(cellI) && visited.get(cellI) == 0)
00563 {
00564
00565 nextCell.append(cellI);
00566
00567
00568
00569
00570
00571
00572 do
00573 {
00574 label currentCell = nextCell.removeHead();
00575
00576 if (visited.get(currentCell) == 0)
00577 {
00578 visited.set(currentCell, 1);
00579
00580
00581 newOrder[cellInOrder] = currentCell;
00582 cellInOrder++;
00583
00584
00585 label i0 = (currentCell == 0 ? 0 : offsets[currentCell-1]);
00586 label i1 = offsets[currentCell];
00587
00588 for (label i = i0; i < i1; i++)
00589 {
00590 label nbr = cellCellAddressing.m()[i];
00591
00592 if (!cellRemoved(nbr) && visited.get(nbr) == 0)
00593 {
00594
00595 nextCell.append(nbr);
00596 }
00597 }
00598 }
00599 }
00600 while (nextCell.size());
00601 }
00602 }
00603
00604
00605
00606 newOrder.setSize(cellInOrder);
00607
00608
00609 oldToNew = invert(cellCellAddressing.size(), newOrder);
00610
00611 return cellInOrder;
00612 }
00613
00614
00615
00616
00617
00618 void Foam::polyTopoChange::getFaceOrder
00619 (
00620 const label nActiveFaces,
00621 const labelList& cellFaces,
00622 const labelList& cellFaceOffsets,
00623
00624 labelList& oldToNew,
00625 labelList& patchSizes,
00626 labelList& patchStarts
00627 ) const
00628 {
00629 oldToNew.setSize(faceOwner_.size());
00630 oldToNew = -1;
00631
00632
00633 label newFaceI = 0;
00634
00635 forAll(cellMap_, cellI)
00636 {
00637 label startOfCell = cellFaceOffsets[cellI];
00638 label nFaces = cellFaceOffsets[cellI+1] - startOfCell;
00639
00640
00641 SortableList<label> nbr(nFaces);
00642
00643 for (label i = 0; i < nFaces; i++)
00644 {
00645 label faceI = cellFaces[startOfCell + i];
00646
00647 label nbrCellI = faceNeighbour_[faceI];
00648
00649 if (faceI >= nActiveFaces)
00650 {
00651
00652 nbr[i] = -1;
00653 }
00654 else if (nbrCellI != -1)
00655 {
00656
00657 if (nbrCellI == cellI)
00658 {
00659 nbrCellI = faceOwner_[faceI];
00660 }
00661
00662 if (cellI < nbrCellI)
00663 {
00664
00665 nbr[i] = nbrCellI;
00666 }
00667 else
00668 {
00669
00670 nbr[i] = -1;
00671 }
00672 }
00673 else
00674 {
00675
00676 nbr[i] = -1;
00677 }
00678 }
00679
00680 nbr.sort();
00681
00682 forAll(nbr, i)
00683 {
00684 if (nbr[i] != -1)
00685 {
00686 oldToNew[cellFaces[startOfCell + nbr.indices()[i]]] =
00687 newFaceI++;
00688 }
00689 }
00690 }
00691
00692
00693
00694 patchStarts.setSize(nPatches_);
00695 patchStarts = 0;
00696 patchSizes.setSize(nPatches_);
00697 patchSizes = 0;
00698
00699 patchStarts[0] = newFaceI;
00700
00701 for (label faceI = 0; faceI < nActiveFaces; faceI++)
00702 {
00703 if (region_[faceI] >= 0)
00704 {
00705 patchSizes[region_[faceI]]++;
00706 }
00707 }
00708
00709 label faceI = patchStarts[0];
00710
00711 forAll(patchStarts, patchI)
00712 {
00713 patchStarts[patchI] = faceI;
00714 faceI += patchSizes[patchI];
00715 }
00716
00717
00718
00719
00720
00721
00722
00723 labelList workPatchStarts(patchStarts);
00724
00725 for (label faceI = 0; faceI < nActiveFaces; faceI++)
00726 {
00727 if (region_[faceI] >= 0)
00728 {
00729 oldToNew[faceI] = workPatchStarts[region_[faceI]]++;
00730 }
00731 }
00732
00733
00734 for (label faceI = nActiveFaces; faceI < oldToNew.size(); faceI++)
00735 {
00736 oldToNew[faceI] = faceI;
00737 }
00738
00739
00740 forAll(oldToNew, faceI)
00741 {
00742 if (oldToNew[faceI] == -1)
00743 {
00744 FatalErrorIn
00745 (
00746 "polyTopoChange::getFaceOrder"
00747 "(const label, const labelList&, const labelList&)"
00748 " const"
00749 ) << "Did not determine new position"
00750 << " for face " << faceI
00751 << abort(FatalError);
00752 }
00753 }
00754 }
00755
00756
00757
00758 void Foam::polyTopoChange::reorderCompactFaces
00759 (
00760 const label newSize,
00761 const labelList& oldToNew
00762 )
00763 {
00764 reorder(oldToNew, faces_);
00765 faces_.setCapacity(newSize);
00766
00767 reorder(oldToNew, region_);
00768 region_.setCapacity(newSize);
00769
00770 reorder(oldToNew, faceOwner_);
00771 faceOwner_.setCapacity(newSize);
00772
00773 reorder(oldToNew, faceNeighbour_);
00774 faceNeighbour_.setCapacity(newSize);
00775
00776
00777 reorder(oldToNew, faceMap_);
00778 faceMap_.setCapacity(newSize);
00779
00780 renumberReverseMap(oldToNew, reverseFaceMap_);
00781
00782 renumberKey(oldToNew, faceFromPoint_);
00783 renumberKey(oldToNew, faceFromEdge_);
00784 inplaceReorder(oldToNew, flipFaceFlux_);
00785 flipFaceFlux_.setCapacity(newSize);
00786
00787 reorder(oldToNew, faceZone_);
00788 faceZone_.setCapacity(newSize);
00789
00790 inplaceReorder(oldToNew, faceZoneFlip_);
00791 faceZoneFlip_.setCapacity(newSize);
00792 }
00793
00794
00795
00796
00797
00798
00799 void Foam::polyTopoChange::compact
00800 (
00801 const bool orderCells,
00802 const bool orderPoints,
00803 label& nInternalPoints,
00804 labelList& patchSizes,
00805 labelList& patchStarts
00806 )
00807 {
00808 points_.shrink();
00809 pointMap_.shrink();
00810 reversePointMap_.shrink();
00811 pointZone_.shrink();
00812
00813 faces_.shrink();
00814 region_.shrink();
00815 faceOwner_.shrink();
00816 faceNeighbour_.shrink();
00817 faceMap_.shrink();
00818 reverseFaceMap_.shrink();
00819 faceZone_.shrink();
00820
00821 cellMap_.shrink();
00822 reverseCellMap_.shrink();
00823 cellZone_.shrink();
00824
00825
00826
00827 label nActivePoints = 0;
00828 {
00829 labelList localPointMap(points_.size(), -1);
00830 label newPointI = 0;
00831
00832 if (!orderPoints)
00833 {
00834 nInternalPoints = -1;
00835
00836 forAll(points_, pointI)
00837 {
00838 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
00839 {
00840 localPointMap[pointI] = newPointI++;
00841 }
00842 }
00843 nActivePoints = newPointI;
00844 }
00845 else
00846 {
00847 forAll(points_, pointI)
00848 {
00849 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
00850 {
00851 nActivePoints++;
00852 }
00853 }
00854
00855
00856 forAll(faceOwner_, faceI)
00857 {
00858 if
00859 (
00860 !faceRemoved(faceI)
00861 && faceOwner_[faceI] >= 0
00862 && faceNeighbour_[faceI] < 0
00863 )
00864 {
00865
00866 const face& f = faces_[faceI];
00867
00868 forAll(f, fp)
00869 {
00870 label pointI = f[fp];
00871
00872 if (localPointMap[pointI] == -1)
00873 {
00874 if
00875 (
00876 pointRemoved(pointI)
00877 || retiredPoints_.found(pointI)
00878 )
00879 {
00880 FatalErrorIn("polyTopoChange::compact(..)")
00881 << "Removed or retired point " << pointI
00882 << " in face " << f
00883 << " at position " << faceI << endl
00884 << "Probably face has not been adapted for"
00885 << " removed points." << abort(FatalError);
00886 }
00887 localPointMap[pointI] = newPointI++;
00888 }
00889 }
00890 }
00891 }
00892
00893 label nBoundaryPoints = newPointI;
00894 nInternalPoints = nActivePoints - nBoundaryPoints;
00895
00896
00897 forAll(localPointMap, pointI)
00898 {
00899 if (localPointMap[pointI] != -1)
00900 {
00901 localPointMap[pointI] += nInternalPoints;
00902 }
00903 }
00904
00905 newPointI = 0;
00906
00907
00908 forAll(faceOwner_, faceI)
00909 {
00910 if
00911 (
00912 !faceRemoved(faceI)
00913 && faceOwner_[faceI] >= 0
00914 && faceNeighbour_[faceI] >= 0
00915 )
00916 {
00917
00918 const face& f = faces_[faceI];
00919
00920 forAll(f, fp)
00921 {
00922 label pointI = f[fp];
00923
00924 if (localPointMap[pointI] == -1)
00925 {
00926 if
00927 (
00928 pointRemoved(pointI)
00929 || retiredPoints_.found(pointI)
00930 )
00931 {
00932 FatalErrorIn("polyTopoChange::compact(..)")
00933 << "Removed or retired point " << pointI
00934 << " in face " << f
00935 << " at position " << faceI << endl
00936 << "Probably face has not been adapted for"
00937 << " removed points." << abort(FatalError);
00938 }
00939 localPointMap[pointI] = newPointI++;
00940 }
00941 }
00942 }
00943 }
00944
00945 if (newPointI != nInternalPoints)
00946 {
00947 FatalErrorIn("polyTopoChange::compact(..)")
00948 << "Problem." << abort(FatalError);
00949 }
00950 newPointI = nActivePoints;
00951 }
00952
00953 forAllConstIter(labelHashSet, retiredPoints_, iter)
00954 {
00955 localPointMap[iter.key()] = newPointI++;
00956 }
00957
00958
00959 if (debug)
00960 {
00961 Pout<< "Points : active:" << nActivePoints
00962 << " removed:" << points_.size()-newPointI << endl;
00963 }
00964
00965 reorder(localPointMap, points_);
00966 points_.setCapacity(newPointI);
00967
00968
00969 reorder(localPointMap, pointMap_);
00970 pointMap_.setCapacity(newPointI);
00971 renumberReverseMap(localPointMap, reversePointMap_);
00972
00973 reorder(localPointMap, pointZone_);
00974 pointZone_.setCapacity(newPointI);
00975 renumber(localPointMap, retiredPoints_);
00976
00977
00978 forAll(faces_, faceI)
00979 {
00980 face& f = faces_[faceI];
00981
00982
00983 renumberCompact(localPointMap, f);
00984
00985 if (!faceRemoved(faceI) && f.size() < 3)
00986 {
00987 FatalErrorIn("polyTopoChange::compact(..)")
00988 << "Created illegal face " << f
00989
00990 << " at position:" << faceI
00991 << " when filtering removed points"
00992 << abort(FatalError);
00993 }
00994 }
00995 }
00996
00997
00998
00999 {
01000 labelList localFaceMap(faces_.size(), -1);
01001 label newFaceI = 0;
01002
01003 forAll(faces_, faceI)
01004 {
01005 if (!faceRemoved(faceI) && faceOwner_[faceI] >= 0)
01006 {
01007 localFaceMap[faceI] = newFaceI++;
01008 }
01009 }
01010 nActiveFaces_ = newFaceI;
01011
01012 forAll(faces_, faceI)
01013 {
01014 if (!faceRemoved(faceI) && faceOwner_[faceI] < 0)
01015 {
01016
01017 localFaceMap[faceI] = newFaceI++;
01018 }
01019 }
01020
01021 if (debug)
01022 {
01023 Pout<< "Faces : active:" << nActiveFaces_
01024 << " removed:" << faces_.size()-newFaceI << endl;
01025 }
01026
01027
01028 reorderCompactFaces(newFaceI, localFaceMap);
01029 }
01030
01031
01032 {
01033 labelList localCellMap;
01034 label newCellI;
01035
01036 if (orderCells)
01037 {
01038
01039 CompactListList<label> cellCells;
01040 makeCellCells(nActiveFaces_, cellCells);
01041
01042
01043 newCellI = getCellOrder(cellCells, localCellMap);
01044 }
01045 else
01046 {
01047
01048 localCellMap.setSize(cellMap_.size());
01049 localCellMap = -1;
01050
01051 newCellI = 0;
01052 forAll(cellMap_, cellI)
01053 {
01054 if (!cellRemoved(cellI))
01055 {
01056 localCellMap[cellI] = newCellI++;
01057 }
01058 }
01059 }
01060
01061 if (debug)
01062 {
01063 Pout<< "Cells : active:" << newCellI
01064 << " removed:" << cellMap_.size()-newCellI << endl;
01065 }
01066
01067
01068 if (orderCells || (newCellI != cellMap_.size()))
01069 {
01070 reorder(localCellMap, cellMap_);
01071 cellMap_.setCapacity(newCellI);
01072 renumberReverseMap(localCellMap, reverseCellMap_);
01073
01074 reorder(localCellMap, cellZone_);
01075 cellZone_.setCapacity(newCellI);
01076
01077 renumberKey(localCellMap, cellFromPoint_);
01078 renumberKey(localCellMap, cellFromEdge_);
01079 renumberKey(localCellMap, cellFromFace_);
01080
01081
01082
01083 forAll(faceOwner_, faceI)
01084 {
01085 label own = faceOwner_[faceI];
01086 label nei = faceNeighbour_[faceI];
01087
01088 if (own >= 0)
01089 {
01090
01091 faceOwner_[faceI] = localCellMap[own];
01092
01093 if (nei >= 0)
01094 {
01095
01096 faceNeighbour_[faceI] = localCellMap[nei];
01097
01098
01099 if
01100 (
01101 faceNeighbour_[faceI] >= 0
01102 && faceNeighbour_[faceI] < faceOwner_[faceI]
01103 )
01104 {
01105 faces_[faceI] = faces_[faceI].reverseFace();
01106 Swap(faceOwner_[faceI], faceNeighbour_[faceI]);
01107 flipFaceFlux_[faceI] =
01108 (
01109 flipFaceFlux_[faceI]
01110 ? 0
01111 : 1
01112 );
01113 faceZoneFlip_[faceI] =
01114 (
01115 faceZoneFlip_[faceI]
01116 ? 0
01117 : 1
01118 );
01119 }
01120 }
01121 }
01122 else if (nei >= 0)
01123 {
01124
01125 faceNeighbour_[faceI] = localCellMap[nei];
01126 }
01127 }
01128 }
01129 }
01130
01131
01132 {
01133
01134 labelList cellFaces;
01135 labelList cellFaceOffsets;
01136 makeCells(nActiveFaces_, cellFaces, cellFaceOffsets);
01137
01138
01139 labelList localFaceMap;
01140 getFaceOrder
01141 (
01142 nActiveFaces_,
01143 cellFaces,
01144 cellFaceOffsets,
01145
01146 localFaceMap,
01147 patchSizes,
01148 patchStarts
01149 );
01150
01151
01152 reorderCompactFaces(localFaceMap.size(), localFaceMap);
01153 }
01154 }
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164 Foam::labelList Foam::polyTopoChange::selectFaces
01165 (
01166 const primitiveMesh& mesh,
01167 const labelList& faceLabels,
01168 const bool internalFacesOnly
01169 )
01170 {
01171 label nFaces = 0;
01172
01173 forAll(faceLabels, i)
01174 {
01175 label faceI = faceLabels[i];
01176
01177 if (internalFacesOnly == mesh.isInternalFace(faceI))
01178 {
01179 nFaces++;
01180 }
01181 }
01182
01183 labelList collectedFaces;
01184
01185 if (nFaces == 0)
01186 {
01187
01188
01189 collectedFaces = faceLabels;
01190 }
01191 else
01192 {
01193 collectedFaces.setSize(nFaces);
01194
01195 nFaces = 0;
01196
01197 forAll(faceLabels, i)
01198 {
01199 label faceI = faceLabels[i];
01200
01201 if (internalFacesOnly == mesh.isInternalFace(faceI))
01202 {
01203 collectedFaces[nFaces++] = faceI;
01204 }
01205 }
01206 }
01207
01208 return collectedFaces;
01209 }
01210
01211
01212
01213
01214 void Foam::polyTopoChange::calcPatchPointMap
01215 (
01216 const List<Map<label> >& oldPatchMeshPointMaps,
01217 const polyBoundaryMesh& boundary,
01218 labelListList& patchPointMap
01219 ) const
01220 {
01221 patchPointMap.setSize(boundary.size());
01222
01223 forAll(boundary, patchI)
01224 {
01225 const labelList& meshPoints = boundary[patchI].meshPoints();
01226
01227 const Map<label>& oldMeshPointMap = oldPatchMeshPointMaps[patchI];
01228
01229 labelList& curPatchPointRnb = patchPointMap[patchI];
01230
01231 curPatchPointRnb.setSize(meshPoints.size());
01232
01233 forAll(meshPoints, i)
01234 {
01235 if (meshPoints[i] < pointMap_.size())
01236 {
01237
01238 Map<label>::const_iterator ozmpmIter = oldMeshPointMap.find
01239 (
01240 pointMap_[meshPoints[i]]
01241 );
01242
01243 if (ozmpmIter != oldMeshPointMap.end())
01244 {
01245 curPatchPointRnb[i] = ozmpmIter();
01246 }
01247 else
01248 {
01249 curPatchPointRnb[i] = -1;
01250 }
01251 }
01252 else
01253 {
01254 curPatchPointRnb[i] = -1;
01255 }
01256 }
01257 }
01258 }
01259
01260
01261 void Foam::polyTopoChange::calcFaceInflationMaps
01262 (
01263 const polyMesh& mesh,
01264 List<objectMap>& facesFromPoints,
01265 List<objectMap>& facesFromEdges,
01266 List<objectMap>& facesFromFaces
01267 ) const
01268 {
01269
01270
01271
01272 facesFromPoints.setSize(faceFromPoint_.size());
01273
01274 if (faceFromPoint_.size())
01275 {
01276 label nFacesFromPoints = 0;
01277
01278
01279 forAllConstIter(Map<label>, faceFromPoint_, iter)
01280 {
01281 label newFaceI = iter.key();
01282
01283 if (region_[newFaceI] == -1)
01284 {
01285
01286 facesFromPoints[nFacesFromPoints++] = objectMap
01287 (
01288 newFaceI,
01289 selectFaces
01290 (
01291 mesh,
01292 mesh.pointFaces()[iter()],
01293 true
01294 )
01295 );
01296 }
01297 else
01298 {
01299
01300 facesFromPoints[nFacesFromPoints++] = objectMap
01301 (
01302 newFaceI,
01303 selectFaces
01304 (
01305 mesh,
01306 mesh.pointFaces()[iter()],
01307 false
01308 )
01309 );
01310 }
01311 }
01312 }
01313
01314
01315
01316
01317
01318 facesFromEdges.setSize(faceFromEdge_.size());
01319
01320 if (faceFromEdge_.size())
01321 {
01322 label nFacesFromEdges = 0;
01323
01324
01325 forAllConstIter(Map<label>, faceFromEdge_, iter)
01326 {
01327 label newFaceI = iter.key();
01328
01329 if (region_[newFaceI] == -1)
01330 {
01331
01332 facesFromEdges[nFacesFromEdges++] = objectMap
01333 (
01334 newFaceI,
01335 selectFaces
01336 (
01337 mesh,
01338 mesh.edgeFaces(iter()),
01339 true
01340 )
01341 );
01342 }
01343 else
01344 {
01345
01346 facesFromEdges[nFacesFromEdges++] = objectMap
01347 (
01348 newFaceI,
01349 selectFaces
01350 (
01351 mesh,
01352 mesh.edgeFaces(iter()),
01353 false
01354 )
01355 );
01356 }
01357 }
01358 }
01359
01360
01361
01362
01363
01364 getMergeSets
01365 (
01366 reverseFaceMap_,
01367 faceMap_,
01368 facesFromFaces
01369 );
01370 }
01371
01372
01373 void Foam::polyTopoChange::calcCellInflationMaps
01374 (
01375 const polyMesh& mesh,
01376 List<objectMap>& cellsFromPoints,
01377 List<objectMap>& cellsFromEdges,
01378 List<objectMap>& cellsFromFaces,
01379 List<objectMap>& cellsFromCells
01380 ) const
01381 {
01382 cellsFromPoints.setSize(cellFromPoint_.size());
01383
01384 if (cellFromPoint_.size())
01385 {
01386 label nCellsFromPoints = 0;
01387
01388
01389 forAllConstIter(Map<label>, cellFromPoint_, iter)
01390 {
01391 cellsFromPoints[nCellsFromPoints++] = objectMap
01392 (
01393 iter.key(),
01394 mesh.pointCells()[iter()]
01395 );
01396 }
01397 }
01398
01399
01400 cellsFromEdges.setSize(cellFromEdge_.size());
01401
01402 if (cellFromEdge_.size())
01403 {
01404 label nCellsFromEdges = 0;
01405
01406
01407 forAllConstIter(Map<label>, cellFromEdge_, iter)
01408 {
01409 cellsFromEdges[nCellsFromEdges++] = objectMap
01410 (
01411 iter.key(),
01412 mesh.edgeCells()[iter()]
01413 );
01414 }
01415 }
01416
01417
01418 cellsFromFaces.setSize(cellFromFace_.size());
01419
01420 if (cellFromFace_.size())
01421 {
01422 label nCellsFromFaces = 0;
01423
01424 labelList twoCells(2);
01425
01426
01427 forAllConstIter(Map<label>, cellFromFace_, iter)
01428 {
01429 label oldFaceI = iter();
01430
01431 if (mesh.isInternalFace(oldFaceI))
01432 {
01433 twoCells[0] = mesh.faceOwner()[oldFaceI];
01434 twoCells[1] = mesh.faceNeighbour()[oldFaceI];
01435 cellsFromFaces[nCellsFromFaces++] = objectMap
01436 (
01437 iter.key(),
01438 twoCells
01439 );
01440 }
01441 else
01442 {
01443 cellsFromFaces[nCellsFromFaces++] = objectMap
01444 (
01445 iter.key(),
01446 labelList(1, mesh.faceOwner()[oldFaceI])
01447 );
01448 }
01449 }
01450 }
01451
01452
01453
01454
01455
01456 getMergeSets
01457 (
01458 reverseCellMap_,
01459 cellMap_,
01460 cellsFromCells
01461 );
01462 }
01463
01464
01465 void Foam::polyTopoChange::resetZones
01466 (
01467 const polyMesh& mesh,
01468 polyMesh& newMesh,
01469 labelListList& pointZoneMap,
01470 labelListList& faceZoneFaceMap,
01471 labelListList& cellZoneMap
01472 ) const
01473 {
01474
01475
01476
01477 pointZoneMap.setSize(mesh.pointZones().size());
01478 {
01479 const pointZoneMesh& pointZones = mesh.pointZones();
01480
01481
01482
01483 labelList nPoints(pointZones.size(), 0);
01484
01485 forAll(pointZone_, pointI)
01486 {
01487 label zoneI = pointZone_[pointI];
01488
01489 if (zoneI >= pointZones.size())
01490 {
01491 FatalErrorIn
01492 (
01493 "resetZones(const polyMesh&, polyMesh&, labelListList&"
01494 "labelListList&, labelListList&)"
01495 ) << "Illegal zoneID " << zoneI << " for point "
01496 << pointI << " coord " << mesh.points()[pointI]
01497 << abort(FatalError);
01498 }
01499
01500 if (zoneI >= 0)
01501 {
01502 nPoints[zoneI]++;
01503 }
01504 }
01505
01506
01507
01508 labelListList addressing(pointZones.size());
01509 forAll(addressing, zoneI)
01510 {
01511 addressing[zoneI].setSize(nPoints[zoneI]);
01512 }
01513 nPoints = 0;
01514
01515 forAll(pointZone_, pointI)
01516 {
01517 label zoneI = pointZone_[pointI];
01518 if (zoneI >= 0)
01519 {
01520 addressing[zoneI][nPoints[zoneI]++] = pointI;
01521 }
01522 }
01523
01524 forAll(addressing, zoneI)
01525 {
01526 stableSort(addressing[zoneI]);
01527 }
01528
01529
01530
01531 forAll(addressing, zoneI)
01532 {
01533 const pointZone& oldZone = pointZones[zoneI];
01534 const labelList& newZoneAddr = addressing[zoneI];
01535
01536 labelList& curPzRnb = pointZoneMap[zoneI];
01537 curPzRnb.setSize(newZoneAddr.size());
01538
01539 forAll(newZoneAddr, i)
01540 {
01541 if (newZoneAddr[i] < pointMap_.size())
01542 {
01543 curPzRnb[i] = oldZone.whichPoint(pointMap_[newZoneAddr[i]]);
01544 }
01545 else
01546 {
01547 curPzRnb[i] = -1;
01548 }
01549 }
01550 }
01551
01552
01553 newMesh.pointZones().clearAddressing();
01554 forAll(newMesh.pointZones(), zoneI)
01555 {
01556 if (debug)
01557 {
01558 Pout<< "pointZone:" << zoneI
01559 << " name:" << newMesh.pointZones()[zoneI].name()
01560 << " size:" << addressing[zoneI].size()
01561 << endl;
01562 }
01563
01564 newMesh.pointZones()[zoneI] = addressing[zoneI];
01565 }
01566 }
01567
01568
01569
01570
01571
01572 faceZoneFaceMap.setSize(mesh.faceZones().size());
01573 {
01574 const faceZoneMesh& faceZones = mesh.faceZones();
01575
01576 labelList nFaces(faceZones.size(), 0);
01577
01578 forAll(faceZone_, faceI)
01579 {
01580 label zoneI = faceZone_[faceI];
01581
01582 if (zoneI >= faceZones.size())
01583 {
01584 FatalErrorIn
01585 (
01586 "resetZones(const polyMesh&, polyMesh&, labelListList&"
01587 "labelListList&, labelListList&)"
01588 ) << "Illegal zoneID " << zoneI << " for face "
01589 << faceI
01590 << abort(FatalError);
01591 }
01592 if (zoneI >= 0)
01593 {
01594 nFaces[zoneI]++;
01595 }
01596 }
01597
01598 labelListList addressing(faceZones.size());
01599 boolListList flipMode(faceZones.size());
01600
01601 forAll(addressing, zoneI)
01602 {
01603 addressing[zoneI].setSize(nFaces[zoneI]);
01604 flipMode[zoneI].setSize(nFaces[zoneI]);
01605 }
01606 nFaces = 0;
01607
01608 forAll(faceZone_, faceI)
01609 {
01610 label zoneI = faceZone_[faceI];
01611
01612 if (zoneI >= 0)
01613 {
01614 label index = nFaces[zoneI]++;
01615 addressing[zoneI][index] = faceI;
01616 flipMode[zoneI][index] = faceZoneFlip_[faceI];
01617 }
01618 }
01619
01620 forAll(addressing, zoneI)
01621 {
01622 labelList newToOld;
01623 sortedOrder(addressing[zoneI], newToOld);
01624 {
01625 labelList newAddressing(addressing[zoneI].size());
01626 forAll(newAddressing, i)
01627 {
01628 newAddressing[i] = addressing[zoneI][newToOld[i]];
01629 }
01630 addressing[zoneI].transfer(newAddressing);
01631 }
01632 {
01633 boolList newFlipMode(flipMode[zoneI].size());
01634 forAll(newFlipMode, i)
01635 {
01636 newFlipMode[i] = flipMode[zoneI][newToOld[i]];
01637 }
01638 flipMode[zoneI].transfer(newFlipMode);
01639 }
01640 }
01641
01642
01643
01644 forAll(addressing, zoneI)
01645 {
01646 const faceZone& oldZone = faceZones[zoneI];
01647 const labelList& newZoneAddr = addressing[zoneI];
01648
01649 labelList& curFzFaceRnb = faceZoneFaceMap[zoneI];
01650
01651 curFzFaceRnb.setSize(newZoneAddr.size());
01652
01653 forAll(newZoneAddr, i)
01654 {
01655 if (newZoneAddr[i] < faceMap_.size())
01656 {
01657 curFzFaceRnb[i] =
01658 oldZone.whichFace(faceMap_[newZoneAddr[i]]);
01659 }
01660 else
01661 {
01662 curFzFaceRnb[i] = -1;
01663 }
01664 }
01665 }
01666
01667
01668
01669 newMesh.faceZones().clearAddressing();
01670 forAll(newMesh.faceZones(), zoneI)
01671 {
01672 if (debug)
01673 {
01674 Pout<< "faceZone:" << zoneI
01675 << " name:" << newMesh.faceZones()[zoneI].name()
01676 << " size:" << addressing[zoneI].size()
01677 << endl;
01678 }
01679
01680 newMesh.faceZones()[zoneI].resetAddressing
01681 (
01682 addressing[zoneI],
01683 flipMode[zoneI]
01684 );
01685 }
01686 }
01687
01688
01689
01690
01691
01692 cellZoneMap.setSize(mesh.cellZones().size());
01693 {
01694 const cellZoneMesh& cellZones = mesh.cellZones();
01695
01696 labelList nCells(cellZones.size(), 0);
01697
01698 forAll(cellZone_, cellI)
01699 {
01700 label zoneI = cellZone_[cellI];
01701
01702 if (zoneI >= cellZones.size())
01703 {
01704 FatalErrorIn
01705 (
01706 "resetZones(const polyMesh&, polyMesh&, labelListList&"
01707 "labelListList&, labelListList&)"
01708 ) << "Illegal zoneID " << zoneI << " for cell "
01709 << cellI << abort(FatalError);
01710 }
01711
01712 if (zoneI >= 0)
01713 {
01714 nCells[zoneI]++;
01715 }
01716 }
01717
01718 labelListList addressing(cellZones.size());
01719 forAll(addressing, zoneI)
01720 {
01721 addressing[zoneI].setSize(nCells[zoneI]);
01722 }
01723 nCells = 0;
01724
01725 forAll(cellZone_, cellI)
01726 {
01727 label zoneI = cellZone_[cellI];
01728
01729 if (zoneI >= 0)
01730 {
01731 addressing[zoneI][nCells[zoneI]++] = cellI;
01732 }
01733 }
01734
01735 forAll(addressing, zoneI)
01736 {
01737 stableSort(addressing[zoneI]);
01738 }
01739
01740
01741
01742 forAll(addressing, zoneI)
01743 {
01744 const cellZone& oldZone = cellZones[zoneI];
01745 const labelList& newZoneAddr = addressing[zoneI];
01746
01747 labelList& curCellRnb = cellZoneMap[zoneI];
01748
01749 curCellRnb.setSize(newZoneAddr.size());
01750
01751 forAll(newZoneAddr, i)
01752 {
01753 if (newZoneAddr[i] < cellMap_.size())
01754 {
01755 curCellRnb[i] =
01756 oldZone.whichCell(cellMap_[newZoneAddr[i]]);
01757 }
01758 else
01759 {
01760 curCellRnb[i] = -1;
01761 }
01762 }
01763 }
01764
01765
01766 newMesh.cellZones().clearAddressing();
01767 forAll(newMesh.cellZones(), zoneI)
01768 {
01769 if (debug)
01770 {
01771 Pout<< "cellZone:" << zoneI
01772 << " name:" << newMesh.cellZones()[zoneI].name()
01773 << " size:" << addressing[zoneI].size()
01774 << endl;
01775 }
01776
01777 newMesh.cellZones()[zoneI] = addressing[zoneI];
01778 }
01779 }
01780 }
01781
01782
01783 void Foam::polyTopoChange::calcFaceZonePointMap
01784 (
01785 const polyMesh& mesh,
01786 const List<Map<label> >& oldFaceZoneMeshPointMaps,
01787 labelListList& faceZonePointMap
01788 ) const
01789 {
01790 const faceZoneMesh& faceZones = mesh.faceZones();
01791
01792 faceZonePointMap.setSize(faceZones.size());
01793
01794 forAll(faceZones, zoneI)
01795 {
01796 const faceZone& newZone = faceZones[zoneI];
01797
01798 const labelList& newZoneMeshPoints = newZone().meshPoints();
01799
01800 const Map<label>& oldZoneMeshPointMap = oldFaceZoneMeshPointMaps[zoneI];
01801
01802 labelList& curFzPointRnb = faceZonePointMap[zoneI];
01803
01804 curFzPointRnb.setSize(newZoneMeshPoints.size());
01805
01806 forAll(newZoneMeshPoints, pointI)
01807 {
01808 if (newZoneMeshPoints[pointI] < pointMap_.size())
01809 {
01810 Map<label>::const_iterator ozmpmIter =
01811 oldZoneMeshPointMap.find
01812 (
01813 pointMap_[newZoneMeshPoints[pointI]]
01814 );
01815
01816 if (ozmpmIter != oldZoneMeshPointMap.end())
01817 {
01818 curFzPointRnb[pointI] = ozmpmIter();
01819 }
01820 else
01821 {
01822 curFzPointRnb[pointI] = -1;
01823 }
01824 }
01825 else
01826 {
01827 curFzPointRnb[pointI] = -1;
01828 }
01829 }
01830 }
01831 }
01832
01833
01834 Foam::face Foam::polyTopoChange::rotateFace
01835 (
01836 const face& f,
01837 const label nPos
01838 )
01839 {
01840 face newF(f.size());
01841
01842 forAll(f, fp)
01843 {
01844 label fp1 = (fp + nPos) % f.size();
01845
01846 if (fp1 < 0)
01847 {
01848 fp1 += f.size();
01849 }
01850
01851 newF[fp1] = f[fp];
01852 }
01853
01854 return newF;
01855 }
01856
01857
01858 void Foam::polyTopoChange::reorderCoupledFaces
01859 (
01860 const bool syncParallel,
01861 const polyBoundaryMesh& boundary,
01862 const labelList& patchStarts,
01863 const labelList& patchSizes,
01864 const pointField& points
01865 )
01866 {
01867
01868
01869 labelList oldToNew(identity(faces_.size()));
01870
01871
01872 labelList rotation(faces_.size(), 0);
01873
01874
01875 forAll(boundary, patchI)
01876 {
01877 if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
01878 {
01879 boundary[patchI].initOrder
01880 (
01881 primitivePatch
01882 (
01883 SubList<face>
01884 (
01885 faces_,
01886 patchSizes[patchI],
01887 patchStarts[patchI]
01888 ),
01889 points
01890 )
01891 );
01892 }
01893 }
01894
01895
01896
01897 bool anyChanged = false;
01898
01899 forAll(boundary, patchI)
01900 {
01901 if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
01902 {
01903 labelList patchFaceMap(patchSizes[patchI], -1);
01904 labelList patchFaceRotation(patchSizes[patchI], 0);
01905
01906 bool changed = boundary[patchI].order
01907 (
01908 primitivePatch
01909 (
01910 SubList<face>
01911 (
01912 faces_,
01913 patchSizes[patchI],
01914 patchStarts[patchI]
01915 ),
01916 points
01917 ),
01918 patchFaceMap,
01919 patchFaceRotation
01920 );
01921
01922 if (changed)
01923 {
01924
01925 label start = patchStarts[patchI];
01926
01927 forAll(patchFaceMap, patchFaceI)
01928 {
01929 oldToNew[patchFaceI + start] =
01930 start + patchFaceMap[patchFaceI];
01931 }
01932
01933 forAll(patchFaceRotation, patchFaceI)
01934 {
01935 rotation[patchFaceI + start] =
01936 patchFaceRotation[patchFaceI];
01937 }
01938
01939 anyChanged = true;
01940 }
01941 }
01942 }
01943
01944 if (syncParallel)
01945 {
01946 reduce(anyChanged, orOp<bool>());
01947 }
01948
01949 if (anyChanged)
01950 {
01951
01952 reorderCompactFaces(oldToNew.size(), oldToNew);
01953
01954
01955 forAll(rotation, faceI)
01956 {
01957 if (rotation[faceI] != 0)
01958 {
01959 faces_[faceI] = rotateFace(faces_[faceI], rotation[faceI]);
01960 }
01961 }
01962 }
01963 }
01964
01965
01966 void Foam::polyTopoChange::compactAndReorder
01967 (
01968 const polyMesh& mesh,
01969 const bool syncParallel,
01970 const bool orderCells,
01971 const bool orderPoints,
01972
01973 label& nInternalPoints,
01974 pointField& newPoints,
01975 labelList& patchSizes,
01976 labelList& patchStarts,
01977 List<objectMap>& pointsFromPoints,
01978 List<objectMap>& facesFromPoints,
01979 List<objectMap>& facesFromEdges,
01980 List<objectMap>& facesFromFaces,
01981 List<objectMap>& cellsFromPoints,
01982 List<objectMap>& cellsFromEdges,
01983 List<objectMap>& cellsFromFaces,
01984 List<objectMap>& cellsFromCells,
01985 List<Map<label> >& oldPatchMeshPointMaps,
01986 labelList& oldPatchNMeshPoints,
01987 labelList& oldPatchStarts,
01988 List<Map<label> >& oldFaceZoneMeshPointMaps
01989 )
01990 {
01991 if (mesh.boundaryMesh().size() != nPatches_)
01992 {
01993 FatalErrorIn("polyTopoChange::compactAndReorder(..)")
01994 << "polyTopoChange was constructed with a mesh with "
01995 << nPatches_ << " patches." << endl
01996 << "The mesh now provided has a different number of patches "
01997 << mesh.boundaryMesh().size()
01998 << " which is illegal" << endl
01999 << abort(FatalError);
02000 }
02001
02002
02003
02004 compact(orderCells, orderPoints, nInternalPoints, patchSizes, patchStarts);
02005
02006
02007
02008 newPoints.transfer(points_);
02009
02010
02011 reorderCoupledFaces
02012 (
02013 syncParallel,
02014 mesh.boundaryMesh(),
02015 patchStarts,
02016 patchSizes,
02017 newPoints
02018 );
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033 getMergeSets
02034 (
02035 reversePointMap_,
02036 pointMap_,
02037 pointsFromPoints
02038 );
02039
02040 calcFaceInflationMaps
02041 (
02042 mesh,
02043 facesFromPoints,
02044 facesFromEdges,
02045 facesFromFaces
02046 );
02047
02048 calcCellInflationMaps
02049 (
02050 mesh,
02051 cellsFromPoints,
02052 cellsFromEdges,
02053 cellsFromFaces,
02054 cellsFromCells
02055 );
02056
02057
02058 {
02059 faceFromPoint_.clearStorage();
02060 faceFromEdge_.clearStorage();
02061
02062 cellFromPoint_.clearStorage();
02063 cellFromEdge_.clearStorage();
02064 cellFromFace_.clearStorage();
02065 }
02066
02067
02068 const polyBoundaryMesh& boundary = mesh.boundaryMesh();
02069
02070
02071 oldPatchMeshPointMaps.setSize(boundary.size());
02072 oldPatchNMeshPoints.setSize(boundary.size());
02073 oldPatchStarts.setSize(boundary.size());
02074
02075 forAll(boundary, patchI)
02076 {
02077
02078 oldPatchMeshPointMaps[patchI] = boundary[patchI].meshPointMap();
02079 oldPatchNMeshPoints[patchI] = boundary[patchI].meshPoints().size();
02080 oldPatchStarts[patchI] = boundary[patchI].start();
02081 }
02082
02083
02084
02085
02086 oldFaceZoneMeshPointMaps.setSize(mesh.faceZones().size());
02087
02088 forAll(mesh.faceZones(), zoneI)
02089 {
02090 const faceZone& oldZone = mesh.faceZones()[zoneI];
02091
02092 oldFaceZoneMeshPointMaps[zoneI] = oldZone().meshPointMap();
02093 }
02094 }
02095
02096
02097
02098
02099
02100 Foam::polyTopoChange::polyTopoChange(const label nPatches, const bool strict)
02101 :
02102 strict_(strict),
02103 nPatches_(nPatches),
02104 points_(0),
02105 pointMap_(0),
02106 reversePointMap_(0),
02107 pointZone_(0),
02108 retiredPoints_(0),
02109 faces_(0),
02110 region_(0),
02111 faceOwner_(0),
02112 faceNeighbour_(0),
02113 faceMap_(0),
02114 reverseFaceMap_(0),
02115 faceFromPoint_(0),
02116 faceFromEdge_(0),
02117 flipFaceFlux_(0),
02118 faceZone_(0),
02119 faceZoneFlip_(0),
02120 nActiveFaces_(0),
02121 cellMap_(0),
02122 reverseCellMap_(0),
02123 cellFromPoint_(0),
02124 cellFromEdge_(0),
02125 cellFromFace_(0),
02126 cellZone_(0)
02127 {}
02128
02129
02130
02131 Foam::polyTopoChange::polyTopoChange
02132 (
02133 const polyMesh& mesh,
02134 const bool strict
02135 )
02136 :
02137 strict_(strict),
02138 nPatches_(0),
02139 points_(0),
02140 pointMap_(0),
02141 reversePointMap_(0),
02142 pointZone_(0),
02143 retiredPoints_(0),
02144 faces_(0),
02145 region_(0),
02146 faceOwner_(0),
02147 faceNeighbour_(0),
02148 faceMap_(0),
02149 reverseFaceMap_(0),
02150 faceFromPoint_(0),
02151 faceFromEdge_(0),
02152 flipFaceFlux_(0),
02153 faceZone_(0),
02154 faceZoneFlip_(0),
02155 nActiveFaces_(0),
02156 cellMap_(0),
02157 reverseCellMap_(0),
02158 cellFromPoint_(0),
02159 cellFromEdge_(0),
02160 cellFromFace_(0),
02161 cellZone_(0)
02162 {
02163 addMesh
02164 (
02165 mesh,
02166 identity(mesh.boundaryMesh().size()),
02167 identity(mesh.pointZones().size()),
02168 identity(mesh.faceZones().size()),
02169 identity(mesh.cellZones().size())
02170 );
02171 }
02172
02173
02174
02175
02176 void Foam::polyTopoChange::clear()
02177 {
02178 points_.clearStorage();
02179 pointMap_.clearStorage();
02180 reversePointMap_.clearStorage();
02181 pointZone_.clearStorage();
02182 retiredPoints_.clearStorage();
02183
02184 faces_.clearStorage();
02185 region_.clearStorage();
02186 faceOwner_.clearStorage();
02187 faceNeighbour_.clearStorage();
02188 faceMap_.clearStorage();
02189 reverseFaceMap_.clearStorage();
02190 faceFromPoint_.clearStorage();
02191 faceFromEdge_.clearStorage();
02192 flipFaceFlux_.clearStorage();
02193 faceZone_.clearStorage();
02194 faceZoneFlip_.clearStorage();
02195 nActiveFaces_ = 0;
02196
02197 cellMap_.clearStorage();
02198 reverseCellMap_.clearStorage();
02199 cellZone_.clearStorage();
02200 cellFromPoint_.clearStorage();
02201 cellFromEdge_.clearStorage();
02202 cellFromFace_.clearStorage();
02203 }
02204
02205
02206 void Foam::polyTopoChange::addMesh
02207 (
02208 const polyMesh& mesh,
02209 const labelList& patchMap,
02210 const labelList& pointZoneMap,
02211 const labelList& faceZoneMap,
02212 const labelList& cellZoneMap
02213 )
02214 {
02215 label maxRegion = nPatches_ - 1;
02216 forAll(patchMap, i)
02217 {
02218 maxRegion = max(maxRegion, patchMap[i]);
02219 }
02220 nPatches_ = maxRegion + 1;
02221
02222
02223
02224 {
02225 const pointField& points = mesh.points();
02226 const pointZoneMesh& pointZones = mesh.pointZones();
02227
02228
02229 points_.setCapacity(points_.size() + points.size());
02230 pointMap_.setCapacity(pointMap_.size() + points.size());
02231 reversePointMap_.setCapacity(reversePointMap_.size() + points.size());
02232 pointZone_.setCapacity(pointZone_.size() + points.size());
02233
02234
02235 labelList newZoneID(points.size(), -1);
02236
02237 forAll(pointZones, zoneI)
02238 {
02239 const labelList& pointLabels = pointZones[zoneI];
02240
02241 forAll(pointLabels, j)
02242 {
02243 newZoneID[pointLabels[j]] = pointZoneMap[zoneI];
02244 }
02245 }
02246
02247
02248 for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
02249 {
02250 addPoint
02251 (
02252 points[pointI],
02253 pointI,
02254 newZoneID[pointI],
02255 true
02256 );
02257 }
02258 }
02259
02260
02261 {
02262 const cellZoneMesh& cellZones = mesh.cellZones();
02263
02264
02265
02266
02267
02268 label nAllCells = mesh.nCells();
02269
02270 cellMap_.setCapacity(cellMap_.size() + nAllCells);
02271 reverseCellMap_.setCapacity(reverseCellMap_.size() + nAllCells);
02272 cellFromPoint_.resize(cellFromPoint_.size() + nAllCells/100);
02273 cellFromEdge_.resize(cellFromEdge_.size() + nAllCells/100);
02274 cellFromFace_.resize(cellFromFace_.size() + nAllCells/100);
02275 cellZone_.setCapacity(cellZone_.size() + nAllCells);
02276
02277
02278
02279 labelList newZoneID(nAllCells, -1);
02280
02281 forAll(cellZones, zoneI)
02282 {
02283 const labelList& cellLabels = cellZones[zoneI];
02284
02285 forAll(cellLabels, j)
02286 {
02287 label cellI = cellLabels[j];
02288
02289 if (newZoneID[cellI] != -1)
02290 {
02291 WarningIn
02292 (
02293 "polyTopoChange::addMesh"
02294 "(const polyMesh&, const labelList&,"
02295 "const labelList&, const labelList&,"
02296 "const labelList&)"
02297 ) << "Cell:" << cellI
02298 << " centre:" << mesh.cellCentres()[cellI]
02299 << " is in two zones:"
02300 << cellZones[newZoneID[cellI]].name()
02301 << " and " << cellZones[zoneI].name() << endl
02302 << " This is not supported."
02303 << " Continuing with first zone only." << endl;
02304 }
02305 else
02306 {
02307 newZoneID[cellI] = cellZoneMap[zoneI];
02308 }
02309 }
02310 }
02311
02312
02313 for (label cellI = 0; cellI < nAllCells; cellI++)
02314 {
02315
02316 addCell(-1, -1, -1, cellI, newZoneID[cellI]);
02317 }
02318 }
02319
02320
02321 {
02322 const polyBoundaryMesh& patches = mesh.boundaryMesh();
02323 const faceList& faces = mesh.faces();
02324 const labelList& faceOwner = mesh.faceOwner();
02325 const labelList& faceNeighbour = mesh.faceNeighbour();
02326 const faceZoneMesh& faceZones = mesh.faceZones();
02327
02328
02329 label nAllFaces = mesh.faces().size();
02330
02331 faces_.setCapacity(faces_.size() + nAllFaces);
02332 region_.setCapacity(region_.size() + nAllFaces);
02333 faceOwner_.setCapacity(faceOwner_.size() + nAllFaces);
02334 faceNeighbour_.setCapacity(faceNeighbour_.size() + nAllFaces);
02335 faceMap_.setCapacity(faceMap_.size() + nAllFaces);
02336 reverseFaceMap_.setCapacity(reverseFaceMap_.size() + nAllFaces);
02337 faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/100);
02338 faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/100);
02339 flipFaceFlux_.setCapacity(flipFaceFlux_.size() + nAllFaces);
02340 faceZone_.setCapacity(faceZone_.size() + nAllFaces);
02341 faceZoneFlip_.setCapacity(faceZoneFlip_.size() + nAllFaces);
02342
02343
02344
02345 labelList newZoneID(nAllFaces, -1);
02346 boolList zoneFlip(nAllFaces, false);
02347
02348 forAll(faceZones, zoneI)
02349 {
02350 const labelList& faceLabels = faceZones[zoneI];
02351 const boolList& flipMap = faceZones[zoneI].flipMap();
02352
02353 forAll(faceLabels, j)
02354 {
02355 newZoneID[faceLabels[j]] = faceZoneMap[zoneI];
02356 zoneFlip[faceLabels[j]] = flipMap[j];
02357 }
02358 }
02359
02360
02361
02362
02363 for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
02364 {
02365 addFace
02366 (
02367 faces[faceI],
02368 faceOwner[faceI],
02369 faceNeighbour[faceI],
02370 -1,
02371 -1,
02372 faceI,
02373 false,
02374 -1,
02375 newZoneID[faceI],
02376 zoneFlip[faceI]
02377 );
02378 }
02379
02380
02381 forAll(patches, patchI)
02382 {
02383 const polyPatch& pp = patches[patchI];
02384
02385 if (pp.start() != faces_.size())
02386 {
02387 FatalErrorIn
02388 (
02389 "polyTopoChange::polyTopoChange"
02390 "(const polyMesh& mesh, const bool strict)"
02391 ) << "Problem : "
02392 << "Patch " << pp.name() << " starts at " << pp.start()
02393 << endl
02394 << "Current face counter at " << faces_.size() << endl
02395 << "Are patches in incremental order?"
02396 << abort(FatalError);
02397 }
02398 forAll(pp, patchFaceI)
02399 {
02400 label faceI = pp.start() + patchFaceI;
02401
02402 addFace
02403 (
02404 faces[faceI],
02405 faceOwner[faceI],
02406 -1,
02407 -1,
02408 -1,
02409 faceI,
02410 false,
02411 patchMap[patchI],
02412 newZoneID[faceI],
02413 zoneFlip[faceI]
02414 );
02415 }
02416 }
02417 }
02418 }
02419
02420
02421 void Foam::polyTopoChange::setCapacity
02422 (
02423 const label nPoints,
02424 const label nFaces,
02425 const label nCells
02426 )
02427 {
02428 points_.setCapacity(nPoints);
02429 pointMap_.setCapacity(nPoints);
02430 reversePointMap_.setCapacity(nPoints);
02431 pointZone_.setCapacity(nPoints);
02432
02433 faces_.setCapacity(nFaces);
02434 region_.setCapacity(nFaces);
02435 faceOwner_.setCapacity(nFaces);
02436 faceNeighbour_.setCapacity(nFaces);
02437 faceMap_.setCapacity(nFaces);
02438 reverseFaceMap_.setCapacity(nFaces);
02439 faceFromPoint_.resize(faceFromPoint_.size() + nFaces/100);
02440 faceFromEdge_.resize(faceFromEdge_.size() + nFaces/100);
02441 flipFaceFlux_.setCapacity(nFaces);
02442 faceZone_.setCapacity(nFaces);
02443 faceZoneFlip_.setCapacity(nFaces);
02444
02445 cellMap_.setCapacity(nCells);
02446 reverseCellMap_.setCapacity(nCells);
02447 cellFromPoint_.resize(cellFromPoint_.size() + nCells/100);
02448 cellFromEdge_.resize(cellFromEdge_.size() + nCells/100);
02449 cellFromFace_.resize(cellFromFace_.size() + nCells/100);
02450 cellZone_.setCapacity(nCells);
02451 }
02452
02453
02454 Foam::label Foam::polyTopoChange::setAction(const topoAction& action)
02455 {
02456 if (isType<polyAddPoint>(action))
02457 {
02458 const polyAddPoint& pap = refCast<const polyAddPoint>(action);
02459
02460 return addPoint
02461 (
02462 pap.newPoint(),
02463 pap.masterPointID(),
02464 pap.zoneID(),
02465 pap.inCell()
02466 );
02467 }
02468 else if (isType<polyModifyPoint>(action))
02469 {
02470 const polyModifyPoint& pmp = refCast<const polyModifyPoint>(action);
02471
02472 modifyPoint
02473 (
02474 pmp.pointID(),
02475 pmp.newPoint(),
02476 pmp.zoneID(),
02477 pmp.inCell()
02478 );
02479
02480 return -1;
02481 }
02482 else if (isType<polyRemovePoint>(action))
02483 {
02484 const polyRemovePoint& prp = refCast<const polyRemovePoint>(action);
02485
02486 removePoint(prp.pointID(), prp.mergePointID());
02487
02488 return -1;
02489 }
02490 else if (isType<polyAddFace>(action))
02491 {
02492 const polyAddFace& paf = refCast<const polyAddFace>(action);
02493
02494 return addFace
02495 (
02496 paf.newFace(),
02497 paf.owner(),
02498 paf.neighbour(),
02499 paf.masterPointID(),
02500 paf.masterEdgeID(),
02501 paf.masterFaceID(),
02502 paf.flipFaceFlux(),
02503 paf.patchID(),
02504 paf.zoneID(),
02505 paf.zoneFlip()
02506 );
02507 }
02508 else if (isType<polyModifyFace>(action))
02509 {
02510 const polyModifyFace& pmf = refCast<const polyModifyFace>(action);
02511
02512 modifyFace
02513 (
02514 pmf.newFace(),
02515 pmf.faceID(),
02516 pmf.owner(),
02517 pmf.neighbour(),
02518 pmf.flipFaceFlux(),
02519 pmf.patchID(),
02520 pmf.zoneID(),
02521 pmf.zoneFlip()
02522 );
02523
02524 return -1;
02525 }
02526 else if (isType<polyRemoveFace>(action))
02527 {
02528 const polyRemoveFace& prf = refCast<const polyRemoveFace>(action);
02529
02530 removeFace(prf.faceID(), prf.mergeFaceID());
02531
02532 return -1;
02533 }
02534 else if (isType<polyAddCell>(action))
02535 {
02536 const polyAddCell& pac = refCast<const polyAddCell>(action);
02537
02538 return addCell
02539 (
02540 pac.masterPointID(),
02541 pac.masterEdgeID(),
02542 pac.masterFaceID(),
02543 pac.masterCellID(),
02544 pac.zoneID()
02545 );
02546 }
02547 else if (isType<polyModifyCell>(action))
02548 {
02549 const polyModifyCell& pmc = refCast<const polyModifyCell>(action);
02550
02551 if (pmc.removeFromZone())
02552 {
02553 modifyCell(pmc.cellID(), -1);
02554 }
02555 else
02556 {
02557 modifyCell(pmc.cellID(), pmc.zoneID());
02558 }
02559
02560 return -1;
02561 }
02562 else if (isType<polyRemoveCell>(action))
02563 {
02564 const polyRemoveCell& prc = refCast<const polyRemoveCell>(action);
02565
02566 removeCell(prc.cellID(), prc.mergeCellID());
02567
02568 return -1;
02569 }
02570 else
02571 {
02572 FatalErrorIn
02573 (
02574 "label polyTopoChange::setAction(const topoAction& action)"
02575 ) << "Unknown type of topoChange: " << action.type()
02576 << abort(FatalError);
02577
02578
02579 return -1;
02580 }
02581 }
02582
02583
02584 Foam::label Foam::polyTopoChange::addPoint
02585 (
02586 const point& pt,
02587 const label masterPointID,
02588 const label zoneID,
02589 const bool inCell
02590 )
02591 {
02592 label pointI = points_.size();
02593
02594 points_.append(pt);
02595 pointMap_.append(masterPointID);
02596 reversePointMap_.append(pointI);
02597 pointZone_.append(zoneID);
02598
02599 if (!inCell)
02600 {
02601 retiredPoints_.insert(pointI);
02602 }
02603
02604 return pointI;
02605 }
02606
02607
02608 void Foam::polyTopoChange::modifyPoint
02609 (
02610 const label pointI,
02611 const point& pt,
02612 const label newZoneID,
02613 const bool inCell
02614 )
02615 {
02616 if (pointI < 0 || pointI >= points_.size())
02617 {
02618 FatalErrorIn
02619 (
02620 "polyTopoChange::modifyPoint(const label, const point&)"
02621 ) << "illegal point label " << pointI << endl
02622 << "Valid point labels are 0 .. " << points_.size()-1
02623 << abort(FatalError);
02624 }
02625 if (pointRemoved(pointI) || pointMap_[pointI] == -1)
02626 {
02627 FatalErrorIn
02628 (
02629 "polyTopoChange::modifyPoint(const label, const point&)"
02630 ) << "point " << pointI << " already marked for removal"
02631 << abort(FatalError);
02632 }
02633 points_[pointI] = pt;
02634 pointZone_[pointI] = newZoneID;
02635
02636 if (inCell)
02637 {
02638 retiredPoints_.erase(pointI);
02639 }
02640 else
02641 {
02642 retiredPoints_.insert(pointI);
02643 }
02644 }
02645
02646
02647 void Foam::polyTopoChange::movePoints(const pointField& newPoints)
02648 {
02649 if (newPoints.size() != points_.size())
02650 {
02651 FatalErrorIn("polyTopoChange::movePoints(const pointField&)")
02652 << "illegal pointField size." << endl
02653 << "Size:" << newPoints.size() << endl
02654 << "Points in mesh:" << points_.size()
02655 << abort(FatalError);
02656 }
02657
02658 forAll(points_, pointI)
02659 {
02660 points_[pointI] = newPoints[pointI];
02661 }
02662 }
02663
02664
02665 void Foam::polyTopoChange::removePoint
02666 (
02667 const label pointI,
02668 const label mergePointI
02669 )
02670 {
02671 if (pointI < 0 || pointI >= points_.size())
02672 {
02673 FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
02674 << "illegal point label " << pointI << endl
02675 << "Valid point labels are 0 .. " << points_.size()-1
02676 << abort(FatalError);
02677 }
02678
02679 if
02680 (
02681 strict_
02682 && (pointRemoved(pointI) || pointMap_[pointI] == -1)
02683 )
02684 {
02685 FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
02686 << "point " << pointI << " already marked for removal" << nl
02687 << "Point:" << points_[pointI] << " pointMap:" << pointMap_[pointI]
02688 << abort(FatalError);
02689 }
02690
02691 if (pointI == mergePointI)
02692 {
02693 FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
02694 << "Cannot remove/merge point " << pointI << " onto itself."
02695 << abort(FatalError);
02696 }
02697
02698 points_[pointI] = greatPoint;
02699 pointMap_[pointI] = -1;
02700 if (mergePointI >= 0)
02701 {
02702 reversePointMap_[pointI] = -mergePointI-2;
02703 }
02704 else
02705 {
02706 reversePointMap_[pointI] = -1;
02707 }
02708 pointZone_[pointI] = -1;
02709 retiredPoints_.erase(pointI);
02710 }
02711
02712
02713 Foam::label Foam::polyTopoChange::addFace
02714 (
02715 const face& f,
02716 const label own,
02717 const label nei,
02718 const label masterPointID,
02719 const label masterEdgeID,
02720 const label masterFaceID,
02721 const bool flipFaceFlux,
02722 const label patchID,
02723 const label zoneID,
02724 const bool zoneFlip
02725 )
02726 {
02727
02728 if (debug)
02729 {
02730 checkFace(f, -1, own, nei, patchID, zoneID);
02731 }
02732
02733 label faceI = faces_.size();
02734
02735 faces_.append(f);
02736 region_.append(patchID);
02737 faceOwner_.append(own);
02738 faceNeighbour_.append(nei);
02739
02740 if (masterPointID >= 0)
02741 {
02742 faceMap_.append(-1);
02743 faceFromPoint_.insert(faceI, masterPointID);
02744 }
02745 else if (masterEdgeID >= 0)
02746 {
02747 faceMap_.append(-1);
02748 faceFromEdge_.insert(faceI, masterEdgeID);
02749 }
02750 else if (masterFaceID >= 0)
02751 {
02752 faceMap_.append(masterFaceID);
02753 }
02754 else
02755 {
02756
02757
02758
02759
02760
02761 faceMap_.append(-1);
02762 }
02763 reverseFaceMap_.append(faceI);
02764
02765 flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
02766 faceZone_.append(zoneID);
02767 faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
02768
02769 return faceI;
02770 }
02771
02772
02773 void Foam::polyTopoChange::modifyFace
02774 (
02775 const face& f,
02776 const label faceI,
02777 const label own,
02778 const label nei,
02779 const bool flipFaceFlux,
02780 const label patchID,
02781 const label zoneID,
02782 const bool zoneFlip
02783 )
02784 {
02785
02786 if (debug)
02787 {
02788 checkFace(f, faceI, own, nei, patchID, zoneID);
02789 }
02790
02791 faces_[faceI] = f;
02792 faceOwner_[faceI] = own;
02793 faceNeighbour_[faceI] = nei;
02794 region_[faceI] = patchID;
02795
02796 flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
02797 faceZone_[faceI] = zoneID;
02798 faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
02799 }
02800
02801
02802 void Foam::polyTopoChange::removeFace(const label faceI, const label mergeFaceI)
02803 {
02804 if (faceI < 0 || faceI >= faces_.size())
02805 {
02806 FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
02807 << "illegal face label " << faceI << endl
02808 << "Valid face labels are 0 .. " << faces_.size()-1
02809 << abort(FatalError);
02810 }
02811
02812 if
02813 (
02814 strict_
02815 && (faceRemoved(faceI) || faceMap_[faceI] == -1)
02816 )
02817 {
02818 FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
02819 << "face " << faceI
02820 << " already marked for removal"
02821 << abort(FatalError);
02822 }
02823
02824 faces_[faceI].setSize(0);
02825 region_[faceI] = -1;
02826 faceOwner_[faceI] = -1;
02827 faceNeighbour_[faceI] = -1;
02828 faceMap_[faceI] = -1;
02829 if (mergeFaceI >= 0)
02830 {
02831 reverseFaceMap_[faceI] = -mergeFaceI-2;
02832 }
02833 else
02834 {
02835 reverseFaceMap_[faceI] = -1;
02836 }
02837 faceFromEdge_.erase(faceI);
02838 faceFromPoint_.erase(faceI);
02839 flipFaceFlux_[faceI] = 0;
02840 faceZone_[faceI] = -1;
02841 faceZoneFlip_[faceI] = 0;
02842 }
02843
02844
02845 Foam::label Foam::polyTopoChange::addCell
02846 (
02847 const label masterPointID,
02848 const label masterEdgeID,
02849 const label masterFaceID,
02850 const label masterCellID,
02851 const label zoneID
02852 )
02853 {
02854 label cellI = cellMap_.size();
02855
02856 if (masterPointID >= 0)
02857 {
02858 cellMap_.append(-1);
02859 cellFromPoint_.insert(cellI, masterPointID);
02860 }
02861 else if (masterEdgeID >= 0)
02862 {
02863 cellMap_.append(-1);
02864 cellFromEdge_.insert(cellI, masterEdgeID);
02865 }
02866 else if (masterFaceID >= 0)
02867 {
02868 cellMap_.append(-1);
02869 cellFromFace_.insert(cellI, masterFaceID);
02870 }
02871 else
02872 {
02873 cellMap_.append(masterCellID);
02874 }
02875 reverseCellMap_.append(cellI);
02876 cellZone_.append(zoneID);
02877
02878 return cellI;
02879 }
02880
02881
02882 void Foam::polyTopoChange::modifyCell
02883 (
02884 const label cellI,
02885 const label zoneID
02886 )
02887 {
02888 cellZone_[cellI] = zoneID;
02889 }
02890
02891
02892 void Foam::polyTopoChange::removeCell(const label cellI, const label mergeCellI)
02893 {
02894 if (cellI < 0 || cellI >= cellMap_.size())
02895 {
02896 FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
02897 << "illegal cell label " << cellI << endl
02898 << "Valid cell labels are 0 .. " << cellMap_.size()-1
02899 << abort(FatalError);
02900 }
02901
02902 if (strict_ && cellMap_[cellI] == -2)
02903 {
02904 FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
02905 << "cell " << cellI
02906 << " already marked for removal"
02907 << abort(FatalError);
02908 }
02909
02910 cellMap_[cellI] = -2;
02911 if (mergeCellI >= 0)
02912 {
02913 reverseCellMap_[cellI] = -mergeCellI-2;
02914 }
02915 else
02916 {
02917 reverseCellMap_[cellI] = -1;
02918 }
02919 cellFromPoint_.erase(cellI);
02920 cellFromEdge_.erase(cellI);
02921 cellFromFace_.erase(cellI);
02922 cellZone_[cellI] = -1;
02923 }
02924
02925
02926 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
02927 (
02928 polyMesh& mesh,
02929 const bool inflate,
02930 const bool syncParallel,
02931 const bool orderCells,
02932 const bool orderPoints
02933 )
02934 {
02935 if (debug)
02936 {
02937 Pout<< "polyTopoChange::changeMesh"
02938 << "(polyMesh&, const bool, const bool, const bool, const bool)"
02939 << endl;
02940 }
02941
02942 if (debug)
02943 {
02944 Pout<< "Old mesh:" << nl;
02945 writeMeshStats(mesh, Pout);
02946 }
02947
02948
02949 pointField newPoints;
02950
02951 label nInternalPoints;
02952
02953 labelList patchSizes;
02954 labelList patchStarts;
02955
02956 List<objectMap> pointsFromPoints;
02957 List<objectMap> facesFromPoints;
02958 List<objectMap> facesFromEdges;
02959 List<objectMap> facesFromFaces;
02960 List<objectMap> cellsFromPoints;
02961 List<objectMap> cellsFromEdges;
02962 List<objectMap> cellsFromFaces;
02963 List<objectMap> cellsFromCells;
02964
02965 List<Map<label> > oldPatchMeshPointMaps;
02966 labelList oldPatchNMeshPoints;
02967 labelList oldPatchStarts;
02968 List<Map<label> > oldFaceZoneMeshPointMaps;
02969
02970
02971 compactAndReorder
02972 (
02973 mesh,
02974 syncParallel,
02975 orderCells,
02976 orderPoints,
02977
02978 nInternalPoints,
02979 newPoints,
02980 patchSizes,
02981 patchStarts,
02982 pointsFromPoints,
02983 facesFromPoints,
02984 facesFromEdges,
02985 facesFromFaces,
02986 cellsFromPoints,
02987 cellsFromEdges,
02988 cellsFromFaces,
02989 cellsFromCells,
02990 oldPatchMeshPointMaps,
02991 oldPatchNMeshPoints,
02992 oldPatchStarts,
02993 oldFaceZoneMeshPointMaps
02994 );
02995
02996 const label nOldPoints(mesh.nPoints());
02997 const label nOldFaces(mesh.nFaces());
02998 const label nOldCells(mesh.nCells());
02999
03000
03001
03002
03003
03004
03005
03006 if (inflate)
03007 {
03008
03009
03010 pointField renumberedMeshPoints(newPoints.size());
03011
03012 forAll(pointMap_, newPointI)
03013 {
03014 label oldPointI = pointMap_[newPointI];
03015
03016 if (oldPointI >= 0)
03017 {
03018 renumberedMeshPoints[newPointI] = mesh.points()[oldPointI];
03019 }
03020 else
03021 {
03022 renumberedMeshPoints[newPointI] = vector::zero;
03023 }
03024 }
03025
03026 mesh.resetPrimitives
03027 (
03028 xferMove(renumberedMeshPoints),
03029 faces_.xfer(),
03030 faceOwner_.xfer(),
03031 faceNeighbour_.xfer(),
03032 patchSizes,
03033 patchStarts,
03034 syncParallel
03035 );
03036
03037 mesh.changing(true);
03038 }
03039 else
03040 {
03041
03042 mesh.resetPrimitives
03043 (
03044 xferMove(newPoints),
03045 faces_.xfer(),
03046 faceOwner_.xfer(),
03047 faceNeighbour_.xfer(),
03048 patchSizes,
03049 patchStarts,
03050 syncParallel
03051 );
03052 mesh.changing(true);
03053 }
03054
03055
03056 {
03057 retiredPoints_.clearStorage();
03058 region_.clearStorage();
03059 }
03060
03061
03062 if (debug)
03063 {
03064
03065 label nAdd, nInflate, nMerge, nRemove;
03066 countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
03067 Pout<< "Points:"
03068 << " added(from point):" << nAdd
03069 << " added(from nothing):" << nInflate
03070 << " merged(into other point):" << nMerge
03071 << " removed:" << nRemove
03072 << nl;
03073
03074 countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
03075 Pout<< "Faces:"
03076 << " added(from face):" << nAdd
03077 << " added(inflated):" << nInflate
03078 << " merged(into other face):" << nMerge
03079 << " removed:" << nRemove
03080 << nl;
03081
03082 countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
03083 Pout<< "Cells:"
03084 << " added(from cell):" << nAdd
03085 << " added(inflated):" << nInflate
03086 << " merged(into other cell):" << nMerge
03087 << " removed:" << nRemove
03088 << nl
03089 << endl;
03090 }
03091
03092 if (debug)
03093 {
03094 Pout<< "New mesh:" << nl;
03095 writeMeshStats(mesh, Pout);
03096 }
03097
03098
03099
03100
03101
03102
03103
03104
03105 labelListList pointZoneMap(mesh.pointZones().size());
03106 labelListList faceZoneFaceMap(mesh.faceZones().size());
03107 labelListList cellZoneMap(mesh.cellZones().size());
03108
03109 resetZones(mesh, mesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
03110
03111
03112 {
03113 pointZone_.clearStorage();
03114 faceZone_.clearStorage();
03115 faceZoneFlip_.clearStorage();
03116 cellZone_.clearStorage();
03117 }
03118
03119
03120
03121
03122
03123 labelListList patchPointMap(mesh.boundaryMesh().size());
03124 calcPatchPointMap
03125 (
03126 oldPatchMeshPointMaps,
03127 mesh.boundaryMesh(),
03128 patchPointMap
03129 );
03130
03131
03132 labelListList faceZonePointMap(mesh.faceZones().size());
03133 calcFaceZonePointMap(mesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
03134
03135 labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
03136
03137 return autoPtr<mapPolyMesh>
03138 (
03139 new mapPolyMesh
03140 (
03141 mesh,
03142 nOldPoints,
03143 nOldFaces,
03144 nOldCells,
03145
03146 pointMap_,
03147 pointsFromPoints,
03148
03149 faceMap_,
03150 facesFromPoints,
03151 facesFromEdges,
03152 facesFromFaces,
03153
03154 cellMap_,
03155 cellsFromPoints,
03156 cellsFromEdges,
03157 cellsFromFaces,
03158 cellsFromCells,
03159
03160 reversePointMap_,
03161 reverseFaceMap_,
03162 reverseCellMap_,
03163
03164 flipFaceFluxSet,
03165
03166 patchPointMap,
03167
03168 pointZoneMap,
03169
03170 faceZonePointMap,
03171 faceZoneFaceMap,
03172 cellZoneMap,
03173
03174 newPoints,
03175 oldPatchStarts,
03176 oldPatchNMeshPoints,
03177 true
03178 )
03179 );
03180
03181
03182
03183 }
03184
03185
03186 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
03187 (
03188 autoPtr<fvMesh>& newMeshPtr,
03189 const IOobject& io,
03190 const fvMesh& mesh,
03191 const bool syncParallel,
03192 const bool orderCells,
03193 const bool orderPoints
03194 )
03195 {
03196 if (debug)
03197 {
03198 Pout<< "polyTopoChange::changeMesh"
03199 << "(autoPtr<fvMesh>&, const IOobject&, const fvMesh&"
03200 << ", const bool, const bool, const bool)"
03201 << endl;
03202 }
03203
03204 if (debug)
03205 {
03206 Pout<< "Old mesh:" << nl;
03207 writeMeshStats(mesh, Pout);
03208 }
03209
03210
03211 pointField newPoints;
03212
03213 label nInternalPoints;
03214
03215 labelList patchSizes;
03216 labelList patchStarts;
03217
03218 List<objectMap> pointsFromPoints;
03219 List<objectMap> facesFromPoints;
03220 List<objectMap> facesFromEdges;
03221 List<objectMap> facesFromFaces;
03222 List<objectMap> cellsFromPoints;
03223 List<objectMap> cellsFromEdges;
03224 List<objectMap> cellsFromFaces;
03225 List<objectMap> cellsFromCells;
03226
03227
03228 List<Map<label> > oldPatchMeshPointMaps;
03229 labelList oldPatchNMeshPoints;
03230 labelList oldPatchStarts;
03231 List<Map<label> > oldFaceZoneMeshPointMaps;
03232
03233
03234 compactAndReorder
03235 (
03236 mesh,
03237 syncParallel,
03238 orderCells,
03239 orderPoints,
03240
03241 nInternalPoints,
03242 newPoints,
03243 patchSizes,
03244 patchStarts,
03245 pointsFromPoints,
03246 facesFromPoints,
03247 facesFromEdges,
03248 facesFromFaces,
03249 cellsFromPoints,
03250 cellsFromEdges,
03251 cellsFromFaces,
03252 cellsFromCells,
03253 oldPatchMeshPointMaps,
03254 oldPatchNMeshPoints,
03255 oldPatchStarts,
03256 oldFaceZoneMeshPointMaps
03257 );
03258
03259 const label nOldPoints(mesh.nPoints());
03260 const label nOldFaces(mesh.nFaces());
03261 const label nOldCells(mesh.nCells());
03262
03263
03264
03265
03266
03267 newMeshPtr.reset
03268 (
03269 new fvMesh
03270 (
03271 io,
03272 xferMove(newPoints),
03273 faces_.xfer(),
03274 faceOwner_.xfer(),
03275 faceNeighbour_.xfer()
03276 )
03277 );
03278 fvMesh& newMesh = newMeshPtr();
03279
03280
03281 {
03282 retiredPoints_.clearStorage();
03283 region_.clearStorage();
03284 }
03285
03286
03287 if (debug)
03288 {
03289
03290 label nAdd, nInflate, nMerge, nRemove;
03291 countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
03292 Pout<< "Points:"
03293 << " added(from point):" << nAdd
03294 << " added(from nothing):" << nInflate
03295 << " merged(into other point):" << nMerge
03296 << " removed:" << nRemove
03297 << nl;
03298
03299 countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
03300 Pout<< "Faces:"
03301 << " added(from face):" << nAdd
03302 << " added(inflated):" << nInflate
03303 << " merged(into other face):" << nMerge
03304 << " removed:" << nRemove
03305 << nl;
03306
03307 countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
03308 Pout<< "Cells:"
03309 << " added(from cell):" << nAdd
03310 << " added(inflated):" << nInflate
03311 << " merged(into other cell):" << nMerge
03312 << " removed:" << nRemove
03313 << nl
03314 << endl;
03315 }
03316
03317
03318 {
03319 const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
03320
03321 List<polyPatch*> newBoundary(oldPatches.size());
03322
03323 forAll(oldPatches, patchI)
03324 {
03325 newBoundary[patchI] = oldPatches[patchI].clone
03326 (
03327 newMesh.boundaryMesh(),
03328 patchI,
03329 patchSizes[patchI],
03330 patchStarts[patchI]
03331 ).ptr();
03332 }
03333 newMesh.addFvPatches(newBoundary);
03334 }
03335
03336
03337
03338
03339
03340
03341 const pointZoneMesh& oldPointZones = mesh.pointZones();
03342 List<pointZone*> pZonePtrs(oldPointZones.size());
03343 {
03344 forAll(oldPointZones, i)
03345 {
03346 pZonePtrs[i] = new pointZone
03347 (
03348 oldPointZones[i].name(),
03349 labelList(0),
03350 i,
03351 newMesh.pointZones()
03352 );
03353 }
03354 }
03355
03356 const faceZoneMesh& oldFaceZones = mesh.faceZones();
03357 List<faceZone*> fZonePtrs(oldFaceZones.size());
03358 {
03359 forAll(oldFaceZones, i)
03360 {
03361 fZonePtrs[i] = new faceZone
03362 (
03363 oldFaceZones[i].name(),
03364 labelList(0),
03365 boolList(0),
03366 i,
03367 newMesh.faceZones()
03368 );
03369 }
03370 }
03371
03372 const cellZoneMesh& oldCellZones = mesh.cellZones();
03373 List<cellZone*> cZonePtrs(oldCellZones.size());
03374 {
03375 forAll(oldCellZones, i)
03376 {
03377 cZonePtrs[i] = new cellZone
03378 (
03379 oldCellZones[i].name(),
03380 labelList(0),
03381 i,
03382 newMesh.cellZones()
03383 );
03384 }
03385 }
03386
03387 newMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
03388
03389
03390
03391
03392 labelListList pointZoneMap(mesh.pointZones().size());
03393 labelListList faceZoneFaceMap(mesh.faceZones().size());
03394 labelListList cellZoneMap(mesh.cellZones().size());
03395
03396 resetZones(mesh, newMesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
03397
03398
03399 {
03400 pointZone_.clearStorage();
03401 faceZone_.clearStorage();
03402 faceZoneFlip_.clearStorage();
03403 cellZone_.clearStorage();
03404 }
03405
03406
03407
03408
03409 labelListList patchPointMap(newMesh.boundaryMesh().size());
03410 calcPatchPointMap
03411 (
03412 oldPatchMeshPointMaps,
03413 newMesh.boundaryMesh(),
03414 patchPointMap
03415 );
03416
03417
03418 labelListList faceZonePointMap(newMesh.faceZones().size());
03419 calcFaceZonePointMap(newMesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
03420
03421 if (debug)
03422 {
03423 Pout<< "New mesh:" << nl;
03424 writeMeshStats(mesh, Pout);
03425 }
03426
03427 labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
03428
03429 return autoPtr<mapPolyMesh>
03430 (
03431 new mapPolyMesh
03432 (
03433 newMesh,
03434 nOldPoints,
03435 nOldFaces,
03436 nOldCells,
03437
03438 pointMap_,
03439 pointsFromPoints,
03440
03441 faceMap_,
03442 facesFromPoints,
03443 facesFromEdges,
03444 facesFromFaces,
03445
03446 cellMap_,
03447 cellsFromPoints,
03448 cellsFromEdges,
03449 cellsFromFaces,
03450 cellsFromCells,
03451
03452 reversePointMap_,
03453 reverseFaceMap_,
03454 reverseCellMap_,
03455
03456 flipFaceFluxSet,
03457
03458 patchPointMap,
03459
03460 pointZoneMap,
03461
03462 faceZonePointMap,
03463 faceZoneFaceMap,
03464 cellZoneMap,
03465
03466 newPoints,
03467 oldPatchStarts,
03468 oldPatchNMeshPoints,
03469 true
03470 )
03471 );
03472
03473
03474
03475 }
03476
03477
03478