FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

polyTopoChange.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
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 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
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 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00059 
00060 // Renumber with special handling for merged items (marked with <-1)
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 // Renumber and remove -1 elements.
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                 // unchanged
00152             }
00153             else
00154             {
00155                 // Added (from another cell v.s. inflated from face/point)
00156                 nAdd++;
00157             }
00158         }
00159         else if (oldCellI == -1)
00160         {
00161             // Created from nothing
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             // unchanged
00178         }
00179         else if (newCellI == -1)
00180         {
00181             // removed
00182             nRemove++;
00183         }
00184         else
00185         {
00186             // merged into -newCellI-2
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     // Per new cell the number of old cells that have been merged into it
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     // From merged cell to set index
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     // Collect cell labels.
00267     // Each objectMap will have
00268     // - index : new mesh cell label
00269     // - masterObjects : list of old cells that have been merged. Element 0
00270     //                   will be the original destination cell label.
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             // oldCellI was merged into mergeCellI
00283 
00284             label setI = cellToMergeSet[mergeCellI];
00285 
00286             objectMap& mergeSet = cellsFromCells[setI];
00287 
00288             if (mergeSet.masterObjects().empty())
00289             {
00290                 // First occurrence of master cell mergeCellI
00291 
00292                 mergeSet.index() = mergeCellI;
00293                 mergeSet.masterObjects().setSize(nMerged[mergeCellI]);
00294 
00295                 // old master label
00296                 mergeSet.masterObjects()[0] = cellMap[mergeCellI];
00297 
00298                 // old slave label
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             // retired face
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     // Faces per cell
00428     labelList nNbrs(cellMap_.size(), 0);
00429 
00430     // 1. Count faces per cell
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     // 2. Calculate offsets
00445 
00446     cellFaceOffsets[0] = 0;
00447     forAll(nNbrs, cellI)
00448     {
00449         cellFaceOffsets[cellI+1] = cellFaceOffsets[cellI] + nNbrs[cellI];
00450     }
00451 
00452     // 3. Fill faces per cell
00453 
00454     // reset the whole list to use as counter
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     // Last offset points to beyond end of cellFaces.
00475     cellFaces.setSize(cellFaceOffsets[cellMap_.size()]);
00476 }
00477 
00478 
00479 // Create cell-cell addressing. Called after compaction (but before ordering)
00480 // of faces
00481 void Foam::polyTopoChange::makeCellCells
00482 (
00483     const label nActiveFaces,
00484     CompactListList<label>& cellCells
00485 ) const
00486 {
00487     // Neighbours per cell
00488     labelList nNbrs(cellMap_.size(), 0);
00489 
00490     // Overall number of cellCells
00491     label nCellCells = 0;
00492 
00493     // 1. Count neighbours (through internal faces) per cell
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     // 2. Calculate offsets
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     // 3. Fill faces per cell
00519 
00520     // reset the whole list to use as counter
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 // Cell ordering (based on bandCompression).
00538 // Handles removed cells. Returns number of remaining cells.
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     // Fifo buffer for string of cells
00550     SLList<label> nextCell;
00551 
00552     // Whether cell has been done already
00553     PackedBoolList visited(cellCellAddressing.size());
00554 
00555     label cellInOrder = 0;
00556 
00557 
00558     // loop over the cells
00559     forAll (visited, cellI)
00560     {
00561         // find the first non-removed cell that has not been visited yet
00562         if (!cellRemoved(cellI) && visited.get(cellI) == 0)
00563         {
00564             // use this cell as a start
00565             nextCell.append(cellI);
00566 
00567             // loop through the nextCell list. Add the first cell into the
00568             // cell order if it has not already been visited and ask for its
00569             // neighbours. If the neighbour in question has not been visited,
00570             // add it to the end of the nextCell list
00571 
00572             do
00573             {
00574                 label currentCell = nextCell.removeHead();
00575 
00576                 if (visited.get(currentCell) == 0)
00577                 {
00578                     visited.set(currentCell, 1);
00579 
00580                     // add into cellOrder
00581                     newOrder[cellInOrder] = currentCell;
00582                     cellInOrder++;
00583 
00584                     // find if the neighbours have been visited
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                             // not visited, add to the list
00595                             nextCell.append(nbr);
00596                         }
00597                     }
00598                 }
00599             }
00600             while (nextCell.size());
00601         }
00602     }
00603 
00604 
00605     // Now we have new-to-old in newOrder.
00606     newOrder.setSize(cellInOrder);
00607 
00608     // Invert to get old-to-new. Make sure removed (i.e. unmapped) cells are -1.
00609     oldToNew = invert(cellCellAddressing.size(), newOrder);
00610 
00611     return cellInOrder;
00612 }
00613 
00614 
00615 // Determine order for faces:
00616 // - upper-triangular order for internal faces
00617 // - external faces after internal faces and in patch order.
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     // First unassigned face
00633     label newFaceI = 0;
00634 
00635     forAll(cellMap_, cellI)
00636     {
00637         label startOfCell = cellFaceOffsets[cellI];
00638         label nFaces = cellFaceOffsets[cellI+1] - startOfCell;
00639 
00640         // Neighbouring cells
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                 // Retired face.
00652                 nbr[i] = -1;
00653             }
00654             else if (nbrCellI != -1)
00655             {
00656                 // Internal face. Get cell on other side.
00657                 if (nbrCellI == cellI)
00658                 {
00659                     nbrCellI = faceOwner_[faceI];
00660                 }
00661 
00662                 if (cellI < nbrCellI)
00663                 {
00664                     // CellI is master
00665                     nbr[i] = nbrCellI;
00666                 }
00667                 else
00668                 {
00669                     // nbrCell is master. Let it handle this face.
00670                     nbr[i] = -1;
00671                 }
00672             }
00673             else
00674             {
00675                 // External face. Do later.
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     // Pick up all patch faces in patch face order.
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     //if (debug)
00718     //{
00719     //    Pout<< "patchSizes:" << patchSizes << nl
00720     //        << "patchStarts:" << patchStarts << endl;
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     // Retired faces.
00734     for (label faceI = nActiveFaces; faceI < oldToNew.size(); faceI++)
00735     {
00736         oldToNew[faceI] = faceI;
00737     }
00738 
00739     // Check done all faces.
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 // Reorder and compact faces according to map.
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     // Update faceMaps.
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 // Compact all and orders points and faces:
00796 // - points into internal followed by external points
00797 // - internalfaces upper-triangular
00798 // - externalfaces after internal ones.
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     // Compact points
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             // Mark boundary points
00856             forAll(faceOwner_, faceI)
00857             {
00858                 if
00859                 (
00860                    !faceRemoved(faceI)
00861                  && faceOwner_[faceI] >= 0
00862                  && faceNeighbour_[faceI] < 0
00863                 )
00864                 {
00865                     // Valid boundary face
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             // Move the boundary addressing up
00897             forAll(localPointMap, pointI)
00898             {
00899                 if (localPointMap[pointI] != -1)
00900                 {
00901                     localPointMap[pointI] += nInternalPoints;
00902                 }
00903             }
00904 
00905             newPointI = 0;
00906 
00907             // Mark internal points
00908             forAll(faceOwner_, faceI)
00909             {
00910                 if
00911                 (
00912                    !faceRemoved(faceI)
00913                  && faceOwner_[faceI] >= 0
00914                  && faceNeighbour_[faceI] >= 0
00915                 )
00916                 {
00917                     // Valid internal face
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         // Update pointMaps
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         // Use map to relabel face vertices
00978         forAll(faces_, faceI)
00979         {
00980             face& f = faces_[faceI];
00981 
00982             //labelList oldF(f);
00983             renumberCompact(localPointMap, f);
00984 
00985             if (!faceRemoved(faceI) && f.size() < 3)
00986             {
00987                 FatalErrorIn("polyTopoChange::compact(..)")
00988                     << "Created illegal face " << f
00989                     //<< " from face " << oldF
00990                     << " at position:" << faceI
00991                     << " when filtering removed points"
00992                     << abort(FatalError);
00993             }
00994         }
00995     }
00996 
00997 
00998     // Compact faces.
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                 // Retired face
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         // Reorder faces.
01028         reorderCompactFaces(newFaceI, localFaceMap);
01029     }
01030 
01031     // Compact cells.
01032     {
01033         labelList localCellMap;
01034         label newCellI;
01035 
01036         if (orderCells)
01037         {
01038             // Construct cellCell addressing
01039             CompactListList<label> cellCells;
01040             makeCellCells(nActiveFaces_, cellCells);
01041 
01042             // Cell ordering (based on bandCompression). Handles removed cells.
01043             newCellI = getCellOrder(cellCells, localCellMap);
01044         }
01045         else
01046         {
01047             // Compact out removed cells
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         // Renumber -if cells reordered or -if cells removed
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             // Renumber owner/neighbour. Take into account if neighbour suddenly
01082             // gets lower cell than owner.
01083             forAll(faceOwner_, faceI)
01084             {
01085                 label own = faceOwner_[faceI];
01086                 label nei = faceNeighbour_[faceI];
01087 
01088                 if (own >= 0)
01089                 {
01090                     // Update owner
01091                     faceOwner_[faceI] = localCellMap[own];
01092 
01093                     if (nei >= 0)
01094                     {
01095                         // Update neighbour.
01096                         faceNeighbour_[faceI] = localCellMap[nei];
01097 
01098                         // Check if face needs reversing.
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                     // Update neighbour.
01125                     faceNeighbour_[faceI] = localCellMap[nei];
01126                 }
01127             }
01128         }
01129     }
01130 
01131     // Reorder faces into upper-triangular and patch ordering
01132     {
01133         // Create cells (packed storage)
01134         labelList cellFaces;
01135         labelList cellFaceOffsets;
01136         makeCells(nActiveFaces_, cellFaces, cellFaceOffsets);
01137 
01138         // Do upper triangular order and patch sorting
01139         labelList localFaceMap;
01140         getFaceOrder
01141         (
01142             nActiveFaces_,
01143             cellFaces,
01144             cellFaceOffsets,
01145 
01146             localFaceMap,
01147             patchSizes,
01148             patchStarts
01149         );
01150 
01151         // Reorder faces.
01152         reorderCompactFaces(localFaceMap.size(), localFaceMap);
01153     }
01154 }
01155 
01156 
01157 // Find faces to interpolate to create value for new face. Only used if
01158 // face was inflated from edge or point. Internal faces should only be
01159 // created from internal faces, external faces only from external faces
01160 // (and ideally the same patch)
01161 // Is bit problematic if there are no faces to select, i.e. in polyDualMesh
01162 // an internal face can be created from a boundary edge with no internal
01163 // faces connected to it.
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         // Did not find any faces of the correct type so just use any old
01188         // face.
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 // Calculate pointMap per patch (so from patch point label to old patch point
01213 // label)
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                 // Check if old point was part of same patch
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     // Faces inflated from points
01270     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
01271 
01272     facesFromPoints.setSize(faceFromPoint_.size());
01273 
01274     if (faceFromPoint_.size())
01275     {
01276         label nFacesFromPoints = 0;
01277 
01278         // Collect all still existing faces connected to this point.
01279         forAllConstIter(Map<label>, faceFromPoint_, iter)
01280         {
01281             label newFaceI = iter.key();
01282 
01283             if (region_[newFaceI] == -1)
01284             {
01285                 // Get internal faces using point on old mesh
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                 // Get patch faces using point on old mesh
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     // Faces inflated from edges
01316     // ~~~~~~~~~~~~~~~~~~~~~~~~~
01317 
01318     facesFromEdges.setSize(faceFromEdge_.size());
01319 
01320     if (faceFromEdge_.size())
01321     {
01322         label nFacesFromEdges = 0;
01323 
01324         // Collect all still existing faces connected to this edge.
01325         forAllConstIter(Map<label>, faceFromEdge_, iter)
01326         {
01327             label newFaceI = iter.key();
01328 
01329             if (region_[newFaceI] == -1)
01330             {
01331                 // Get internal faces using edge on old mesh
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                 // Get patch faces using edge on old mesh
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     // Faces from face merging
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         // Collect all still existing faces connected to this point.
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         // Collect all still existing faces connected to this point.
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         // Collect all still existing faces connected to this point.
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     // Cells from cell merging
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     // pointZones
01475     // ~~~~~~~~~~
01476 
01477     pointZoneMap.setSize(mesh.pointZones().size());
01478     {
01479         const pointZoneMesh& pointZones = mesh.pointZones();
01480 
01481         // Count points per zone
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         // Distribute points per zone
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         // Sort the addressing
01524         forAll(addressing, zoneI)
01525         {
01526             stableSort(addressing[zoneI]);
01527         }
01528 
01529         // So now we both have old zones and the new addressing.
01530         // Invert the addressing to get pointZoneMap.
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         // Reset the addresing on the zone
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     // faceZones
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         // Sort the addressing
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         // So now we both have old zones and the new addressing.
01643         // Invert the addressing to get faceZoneFaceMap.
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         // Reset the addresing on the zone
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     // cellZones
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         // Sort the addressing
01735         forAll(addressing, zoneI)
01736         {
01737             stableSort(addressing[zoneI]);
01738         }
01739 
01740         // So now we both have old zones and the new addressing.
01741         // Invert the addressing to get cellZoneMap.
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         // Reset the addresing on the zone
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     // Mapping for faces (old to new). Extends over all mesh faces for
01868     // convenience (could be just the external faces)
01869     labelList oldToNew(identity(faces_.size()));
01870 
01871     // Rotation on new faces.
01872     labelList rotation(faces_.size(), 0);
01873 
01874     // Send ordering
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     // Receive and calculate ordering
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                 // Merge patch face reordering into mesh face reordering table
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         // Reorder faces according to oldToNew.
01952         reorderCompactFaces(oldToNew.size(), oldToNew);
01953 
01954         // Rotate faces (rotation is already in new face indices).
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     // Remove any holes from points/faces/cells and sort faces.
02003     // Sets nActiveFaces_.
02004     compact(orderCells, orderPoints, nInternalPoints, patchSizes, patchStarts);
02005 
02006     // Transfer points to pointField. points_ are now cleared!
02007     // Only done since e.g. reorderCoupledFaces requires pointField.
02008     newPoints.transfer(points_);
02009 
02010     // Reorder any coupled faces
02011     reorderCoupledFaces
02012     (
02013         syncParallel,
02014         mesh.boundaryMesh(),
02015         patchStarts,
02016         patchSizes,
02017         newPoints
02018     );
02019 
02020 
02021     // Calculate inflation/merging maps
02022     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
02023     // These are for the new face(/point/cell) the old faces whose value
02024     // needs to be
02025     // averaged/summed to get the new value. These old faces come either from
02026     // merged old faces (face remove into other face),
02027     // the old edgeFaces (inflate from edge) or the old pointFaces (inflate
02028     // from point). As an additional complexity will use only internal faces
02029     // to create new value for internal face and vice versa only patch
02030     // faces to to create patch face value.
02031 
02032     // For point only point merging
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     // Clear inflation info
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     // Grab patch mesh point maps
02071     oldPatchMeshPointMaps.setSize(boundary.size());
02072     oldPatchNMeshPoints.setSize(boundary.size());
02073     oldPatchStarts.setSize(boundary.size());
02074 
02075     forAll(boundary, patchI)
02076     {
02077         // Copy old face zone mesh point maps
02078         oldPatchMeshPointMaps[patchI] = boundary[patchI].meshPointMap();
02079         oldPatchNMeshPoints[patchI] = boundary[patchI].meshPoints().size();
02080         oldPatchStarts[patchI] = boundary[patchI].start();
02081     }
02082 
02083     // Grab old face zone mesh point maps.
02084     // These need to be saved before resetting the mesh and are used
02085     // later on to calculate the faceZone pointMaps.
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 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
02098 
02099 // Construct from components
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 // Construct from components
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 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
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     // Add points
02224     {
02225         const pointField& points = mesh.points();
02226         const pointZoneMesh& pointZones = mesh.pointZones();
02227 
02228         // Extend
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         // Precalc offset zones
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         // Add points in mesh order
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     // Add cells
02261     {
02262         const cellZoneMesh& cellZones = mesh.cellZones();
02263 
02264         // Resize
02265 
02266         // Note: polyMesh does not allow retired cells anymore. So allCells
02267         // always equals nCells
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         // Precalc offset zones
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         // Add cells in mesh order
02313         for (label cellI = 0; cellI < nAllCells; cellI++)
02314         {
02315             // Add cell from cell
02316             addCell(-1, -1, -1, cellI, newZoneID[cellI]);
02317         }
02318     }
02319 
02320     // Add faces
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         // Resize
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         // Precalc offset zones
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         // Add faces in mesh order
02361 
02362         // 1. Internal faces
02363         for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
02364         {
02365             addFace
02366             (
02367                 faces[faceI],
02368                 faceOwner[faceI],
02369                 faceNeighbour[faceI],
02370                 -1,                         // masterPointID
02371                 -1,                         // masterEdgeID
02372                 faceI,                      // masterFaceID
02373                 false,                      // flipFaceFlux
02374                 -1,                         // patchID
02375                 newZoneID[faceI],           // zoneID
02376                 zoneFlip[faceI]             // zoneFlip
02377             );
02378         }
02379 
02380         // 2. Patch faces
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,                         // neighbour
02407                     -1,                         // masterPointID
02408                     -1,                         // masterEdgeID
02409                     faceI,                      // masterFaceID
02410                     false,                      // flipFaceFlux
02411                     patchMap[patchI],           // patchID
02412                     newZoneID[faceI],           // zoneID
02413                     zoneFlip[faceI]             // zoneFlip
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         // Dummy return to keep compiler happy
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     // Check validity
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         // Allow inflate-from-nothing?
02757         //FatalErrorIn("polyTopoChange::addFace")
02758         //    << "Need to specify a master point, edge or face"
02759         //    << "face:" << f << " own:" << own << " nei:" << nei
02760         //    << abort(FatalError);
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     // Check validity
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     // new mesh points
02949     pointField newPoints;
02950     // number of internal points
02951     label nInternalPoints;
02952     // patch slicing
02953     labelList patchSizes;
02954     labelList patchStarts;
02955     // inflate maps
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     // old mesh info
02965     List<Map<label> > oldPatchMeshPointMaps;
02966     labelList oldPatchNMeshPoints;
02967     labelList oldPatchStarts;
02968     List<Map<label> > oldFaceZoneMeshPointMaps;
02969 
02970     // Compact, reorder patch faces and calculate mesh/patch maps.
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     // Change the mesh
03002     // ~~~~~~~~~~~~~~~
03003     // This will invalidate any addressing so better make sure you have
03004     // all the information you need!!!
03005 
03006     if (inflate)
03007     {
03008         // Keep (renumbered) mesh points, store new points in map for inflation
03009         // (appended points (i.e. from nowhere) get value zero)
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         // Set new points.
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     // Clear out primitives
03056     {
03057         retiredPoints_.clearStorage();
03058         region_.clearStorage();
03059     }
03060 
03061 
03062     if (debug)
03063     {
03064         // Some stats on changes
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     // Zones
03100     // ~~~~~
03101 
03102     // Inverse of point/face/cell zone addressing.
03103     // For every preserved point/face/cells in zone give the old position.
03104     // For added points, the index is set to -1
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     // Clear zone info
03112     {
03113         pointZone_.clearStorage();
03114         faceZone_.clearStorage();
03115         faceZoneFlip_.clearStorage();
03116         cellZone_.clearStorage();
03117     }
03118 
03119 
03120     // Patch point renumbering
03121     // For every preserved point on a patch give the old position.
03122     // For added points, the index is set to -1
03123     labelListList patchPointMap(mesh.boundaryMesh().size());
03124     calcPatchPointMap
03125     (
03126         oldPatchMeshPointMaps,
03127         mesh.boundaryMesh(),
03128         patchPointMap
03129     );
03130 
03131     // Create the face zone mesh point renumbering
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,          // if empty signals no inflation.
03175             oldPatchStarts,
03176             oldPatchNMeshPoints,
03177             true                // steal storage.
03178         )
03179     );
03180 
03181     // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
03182     // be invalid.
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     // new mesh points
03211     pointField newPoints;
03212     // number of internal points
03213     label nInternalPoints;
03214     // patch slicing
03215     labelList patchSizes;
03216     labelList patchStarts;
03217     // inflate maps
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     // old mesh info
03228     List<Map<label> > oldPatchMeshPointMaps;
03229     labelList oldPatchNMeshPoints;
03230     labelList oldPatchStarts;
03231     List<Map<label> > oldFaceZoneMeshPointMaps;
03232 
03233     // Compact, reorder patch faces and calculate mesh/patch maps.
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     // Create the mesh
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     // Clear out primitives
03281     {
03282         retiredPoints_.clearStorage();
03283         region_.clearStorage();
03284     }
03285 
03286 
03287     if (debug)
03288     {
03289         // Some stats on changes
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     // Zones
03338     // ~~~~~
03339 
03340     // Start off from empty zones.
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     // Inverse of point/face/cell zone addressing.
03390     // For every preserved point/face/cells in zone give the old position.
03391     // For added points, the index is set to -1
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     // Clear zone info
03399     {
03400         pointZone_.clearStorage();
03401         faceZone_.clearStorage();
03402         faceZoneFlip_.clearStorage();
03403         cellZone_.clearStorage();
03404     }
03405 
03406     // Patch point renumbering
03407     // For every preserved point on a patch give the old position.
03408     // For added points, the index is set to -1
03409     labelListList patchPointMap(newMesh.boundaryMesh().size());
03410     calcPatchPointMap
03411     (
03412         oldPatchMeshPointMaps,
03413         newMesh.boundaryMesh(),
03414         patchPointMap
03415     );
03416 
03417     // Create the face zone mesh point renumbering
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,          // if empty signals no inflation.
03467             oldPatchStarts,
03468             oldPatchNMeshPoints,
03469             true                // steal storage.
03470         )
03471     );
03472 
03473     // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
03474     // be invalid.
03475 }
03476 
03477 
03478 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines