FreeFOAM The Cross-Platform CFD Toolkit
polyTopoChange.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "polyTopoChange.H"
27 #include <OpenFOAM/SortableList.H>
28 #include <OpenFOAM/polyMesh.H>
38 #include <OpenFOAM/objectMap.H>
40 #include <finiteVolume/fvMesh.H>
41 
42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43 
44 namespace Foam
45 {
46  defineTypeNameAndDebug(polyTopoChange, 0);
47 }
48 
49 
50 const Foam::point Foam::polyTopoChange::greatPoint
51 (
52  GREAT,
53  GREAT,
54  GREAT
55 );
56 
57 
58 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
59 
60 // Renumber with special handling for merged items (marked with <-1)
61 void Foam::polyTopoChange::renumberReverseMap
62 (
63  const labelList& map,
64  DynamicList<label>& elems
65 )
66 {
67  forAll(elems, elemI)
68  {
69  label val = elems[elemI];
70 
71  if (val >= 0)
72  {
73  elems[elemI] = map[val];
74  }
75  else if (val < -1)
76  {
77  label mergedVal = -val-2;
78  elems[elemI] = -map[mergedVal]-2;
79  }
80  }
81 }
82 
83 
85 (
86  const labelList& map,
87  labelHashSet& elems
88 )
89 {
90  labelHashSet newElems(elems.size());
91 
92  forAllConstIter(labelHashSet, elems, iter)
93  {
94  label newElem = map[iter.key()];
95 
96  if (newElem >= 0)
97  {
98  newElems.insert(newElem);
99  }
100  }
101 
102  elems.transfer(newElems);
103 }
104 
105 
106 // Renumber and remove -1 elements.
107 void Foam::polyTopoChange::renumberCompact
108 (
109  const labelList& map,
110  labelList& elems
111 )
112 {
113  label newElemI = 0;
114 
115  forAll(elems, elemI)
116  {
117  label newVal = map[elems[elemI]];
118 
119  if (newVal != -1)
120  {
121  elems[newElemI++] = newVal;
122  }
123  }
124  elems.setSize(newElemI);
125 }
126 
127 
128 void Foam::polyTopoChange::countMap
129 (
130  const labelList& map,
131  const labelList& reverseMap,
132  label& nAdd,
133  label& nInflate,
134  label& nMerge,
135  label& nRemove
136 )
137 {
138  nAdd = 0;
139  nInflate = 0;
140  nMerge = 0;
141  nRemove = 0;
142 
143  forAll(map, newCellI)
144  {
145  label oldCellI = map[newCellI];
146 
147  if (oldCellI >= 0)
148  {
149  if (reverseMap[oldCellI] == newCellI)
150  {
151  // unchanged
152  }
153  else
154  {
155  // Added (from another cell v.s. inflated from face/point)
156  nAdd++;
157  }
158  }
159  else if (oldCellI == -1)
160  {
161  // Created from nothing
162  nInflate++;
163  }
164  else
165  {
166  FatalErrorIn("countMap") << "old:" << oldCellI
167  << " new:" << newCellI << abort(FatalError);
168  }
169  }
170 
171  forAll(reverseMap, oldCellI)
172  {
173  label newCellI = reverseMap[oldCellI];
174 
175  if (newCellI >= 0)
176  {
177  // unchanged
178  }
179  else if (newCellI == -1)
180  {
181  // removed
182  nRemove++;
183  }
184  else
185  {
186  // merged into -newCellI-2
187  nMerge++;
188  }
189  }
190 }
191 
192 
193 Foam::labelHashSet Foam::polyTopoChange::getSetIndices
194 (
195  const PackedBoolList& lst
196 )
197 {
198  labelHashSet values(lst.count());
199  forAll(lst, i)
200  {
201  if (lst[i])
202  {
203  values.insert(i);
204  }
205  }
206  return values;
207 }
208 
209 
210 void Foam::polyTopoChange::writeMeshStats(const polyMesh& mesh, Ostream& os)
211 {
212  const polyBoundaryMesh& patches = mesh.boundaryMesh();
213 
214  labelList patchSizes(patches.size());
215  labelList patchStarts(patches.size());
216  forAll(patches, patchI)
217  {
218  patchSizes[patchI] = patches[patchI].size();
219  patchStarts[patchI] = patches[patchI].start();
220  }
221 
222  os << " Points : " << mesh.nPoints() << nl
223  << " Faces : " << mesh.nFaces() << nl
224  << " Cells : " << mesh.nCells() << nl
225  << " PatchSizes : " << patchSizes << nl
226  << " PatchStarts : " << patchStarts << nl
227  << endl;
228 }
229 
230 
231 void Foam::polyTopoChange::getMergeSets
232 (
233  const labelList& reverseCellMap,
234  const labelList& cellMap,
235  List<objectMap>& cellsFromCells
236 )
237 {
238  // Per new cell the number of old cells that have been merged into it
239  labelList nMerged(cellMap.size(), 1);
240 
241  forAll(reverseCellMap, oldCellI)
242  {
243  label newCellI = reverseCellMap[oldCellI];
244 
245  if (newCellI < -1)
246  {
247  label mergeCellI = -newCellI-2;
248 
249  nMerged[mergeCellI]++;
250  }
251  }
252 
253  // From merged cell to set index
254  labelList cellToMergeSet(cellMap.size(), -1);
255 
256  label nSets = 0;
257 
258  forAll(nMerged, cellI)
259  {
260  if (nMerged[cellI] > 1)
261  {
262  cellToMergeSet[cellI] = nSets++;
263  }
264  }
265 
266  // Collect cell labels.
267  // Each objectMap will have
268  // - index : new mesh cell label
269  // - masterObjects : list of old cells that have been merged. Element 0
270  // will be the original destination cell label.
271 
272  cellsFromCells.setSize(nSets);
273 
274  forAll(reverseCellMap, oldCellI)
275  {
276  label newCellI = reverseCellMap[oldCellI];
277 
278  if (newCellI < -1)
279  {
280  label mergeCellI = -newCellI-2;
281 
282  // oldCellI was merged into mergeCellI
283 
284  label setI = cellToMergeSet[mergeCellI];
285 
286  objectMap& mergeSet = cellsFromCells[setI];
287 
288  if (mergeSet.masterObjects().empty())
289  {
290  // First occurrence of master cell mergeCellI
291 
292  mergeSet.index() = mergeCellI;
293  mergeSet.masterObjects().setSize(nMerged[mergeCellI]);
294 
295  // old master label
296  mergeSet.masterObjects()[0] = cellMap[mergeCellI];
297 
298  // old slave label
299  mergeSet.masterObjects()[1] = oldCellI;
300 
301  nMerged[mergeCellI] = 2;
302  }
303  else
304  {
305  mergeSet.masterObjects()[nMerged[mergeCellI]++] = oldCellI;
306  }
307  }
308  }
309 }
310 
311 
312 void Foam::polyTopoChange::checkFace
313 (
314  const face& f,
315  const label faceI,
316  const label own,
317  const label nei,
318  const label patchI,
319  const label zoneI
320 ) const
321 {
322  if (nei == -1)
323  {
324  if (own == -1 && zoneI != -1)
325  {
326  // retired face
327  }
328  else if (patchI == -1 || patchI >= nPatches_)
329  {
331  (
332  "polyTopoChange::checkFace(const face&, const label"
333  ", const label, const label, const label)"
334  ) << "Face has no neighbour (so external) but does not have"
335  << " a valid patch" << nl
336  << "f:" << f
337  << " faceI(-1 if added face):" << faceI
338  << " own:" << own << " nei:" << nei
339  << " patchI:" << patchI << abort(FatalError);
340  }
341  }
342  else
343  {
344  if (patchI != -1)
345  {
347  (
348  "polyTopoChange::checkFace(const face&, const label"
349  ", const label, const label, const label)"
350  ) << "Cannot both have valid patchI and neighbour" << nl
351  << "f:" << f
352  << " faceI(-1 if added face):" << faceI
353  << " own:" << own << " nei:" << nei
354  << " patchI:" << patchI << abort(FatalError);
355  }
356 
357  if (nei <= own)
358  {
360  (
361  "polyTopoChange::checkFace(const face&, const label"
362  ", const label, const label, const label)"
363  ) << "Owner cell label should be less than neighbour cell label"
364  << nl
365  << "f:" << f
366  << " faceI(-1 if added face):" << faceI
367  << " own:" << own << " nei:" << nei
368  << " patchI:" << patchI << abort(FatalError);
369  }
370  }
371 
372  if (f.size() < 3 || findIndex(f, -1) != -1)
373  {
375  (
376  "polyTopoChange::checkFace(const face&, const label"
377  ", const label, const label, const label)"
378  ) << "Illegal vertices in face"
379  << nl
380  << "f:" << f
381  << " faceI(-1 if added face):" << faceI
382  << " own:" << own << " nei:" << nei
383  << " patchI:" << patchI << abort(FatalError);
384  }
385  if (faceI >= 0 && faceI < faces_.size() && faceRemoved(faceI))
386  {
388  (
389  "polyTopoChange::checkFace(const face&, const label"
390  ", const label, const label, const label)"
391  ) << "Face already marked for removal"
392  << nl
393  << "f:" << f
394  << " faceI(-1 if added face):" << faceI
395  << " own:" << own << " nei:" << nei
396  << " patchI:" << patchI << abort(FatalError);
397  }
398  forAll(f, fp)
399  {
400  if (f[fp] < points_.size() && pointRemoved(f[fp]))
401  {
403  (
404  "polyTopoChange::checkFace(const face&, const label"
405  ", const label, const label, const label)"
406  ) << "Face uses removed vertices"
407  << nl
408  << "f:" << f
409  << " faceI(-1 if added face):" << faceI
410  << " own:" << own << " nei:" << nei
411  << " patchI:" << patchI << abort(FatalError);
412  }
413  }
414 }
415 
416 
417 void Foam::polyTopoChange::makeCells
418 (
419  const label nActiveFaces,
420  labelList& cellFaces,
421  labelList& cellFaceOffsets
422 ) const
423 {
424  cellFaces.setSize(2*nActiveFaces);
425  cellFaceOffsets.setSize(cellMap_.size() + 1);
426 
427  // Faces per cell
428  labelList nNbrs(cellMap_.size(), 0);
429 
430  // 1. Count faces per cell
431 
432  for (label faceI = 0; faceI < nActiveFaces; faceI++)
433  {
434  nNbrs[faceOwner_[faceI]]++;
435  }
436  for (label faceI = 0; faceI < nActiveFaces; faceI++)
437  {
438  if (faceNeighbour_[faceI] >= 0)
439  {
440  nNbrs[faceNeighbour_[faceI]]++;
441  }
442  }
443 
444  // 2. Calculate offsets
445 
446  cellFaceOffsets[0] = 0;
447  forAll(nNbrs, cellI)
448  {
449  cellFaceOffsets[cellI+1] = cellFaceOffsets[cellI] + nNbrs[cellI];
450  }
451 
452  // 3. Fill faces per cell
453 
454  // reset the whole list to use as counter
455  nNbrs = 0;
456 
457  for (label faceI = 0; faceI < nActiveFaces; faceI++)
458  {
459  label cellI = faceOwner_[faceI];
460 
461  cellFaces[cellFaceOffsets[cellI] + nNbrs[cellI]++] = faceI;
462  }
463 
464  for (label faceI = 0; faceI < nActiveFaces; faceI++)
465  {
466  label cellI = faceNeighbour_[faceI];
467 
468  if (cellI >= 0)
469  {
470  cellFaces[cellFaceOffsets[cellI] + nNbrs[cellI]++] = faceI;
471  }
472  }
473 
474  // Last offset points to beyond end of cellFaces.
475  cellFaces.setSize(cellFaceOffsets[cellMap_.size()]);
476 }
477 
478 
479 // Create cell-cell addressing. Called after compaction (but before ordering)
480 // of faces
481 void Foam::polyTopoChange::makeCellCells
482 (
483  const label nActiveFaces,
484  CompactListList<label>& cellCells
485 ) const
486 {
487  // Neighbours per cell
488  labelList nNbrs(cellMap_.size(), 0);
489 
490  // Overall number of cellCells
491  label nCellCells = 0;
492 
493  // 1. Count neighbours (through internal faces) per cell
494 
495  for (label faceI = 0; faceI < nActiveFaces; faceI++)
496  {
497  if (faceNeighbour_[faceI] >= 0)
498  {
499  nNbrs[faceOwner_[faceI]]++;
500  nNbrs[faceNeighbour_[faceI]]++;
501  nCellCells += 2;
502  }
503  }
504 
505  cellCells.setSize(cellMap_.size(), nCellCells);
506 
507  // 2. Calculate offsets
508 
509  labelList& offsets = cellCells.offsets();
510 
511  label sumSize = 0;
512  forAll(nNbrs, cellI)
513  {
514  sumSize += nNbrs[cellI];
515  offsets[cellI] = sumSize;
516  }
517 
518  // 3. Fill faces per cell
519 
520  // reset the whole list to use as counter
521  nNbrs = 0;
522 
523  for (label faceI = 0; faceI < nActiveFaces; faceI++)
524  {
525  label nei = faceNeighbour_[faceI];
526 
527  if (nei >= 0)
528  {
529  label own = faceOwner_[faceI];
530  cellCells.m()[cellCells.index(own, nNbrs[own]++)] = nei;
531  cellCells.m()[cellCells.index(nei, nNbrs[nei]++)] = own;
532  }
533  }
534 }
535 
536 
537 // Cell ordering (based on bandCompression).
538 // Handles removed cells. Returns number of remaining cells.
539 Foam::label Foam::polyTopoChange::getCellOrder
540 (
541  const CompactListList<label>& cellCellAddressing,
542  labelList& oldToNew
543 ) const
544 {
545  const labelList& offsets = cellCellAddressing.offsets();
546 
547  labelList newOrder(cellCellAddressing.size());
548 
549  // Fifo buffer for string of cells
550  SLList<label> nextCell;
551 
552  // Whether cell has been done already
553  PackedBoolList visited(cellCellAddressing.size());
554 
555  label cellInOrder = 0;
556 
557 
558  // loop over the cells
559  forAll (visited, cellI)
560  {
561  // find the first non-removed cell that has not been visited yet
562  if (!cellRemoved(cellI) && visited.get(cellI) == 0)
563  {
564  // use this cell as a start
565  nextCell.append(cellI);
566 
567  // loop through the nextCell list. Add the first cell into the
568  // cell order if it has not already been visited and ask for its
569  // neighbours. If the neighbour in question has not been visited,
570  // add it to the end of the nextCell list
571 
572  do
573  {
574  label currentCell = nextCell.removeHead();
575 
576  if (visited.get(currentCell) == 0)
577  {
578  visited.set(currentCell, 1);
579 
580  // add into cellOrder
581  newOrder[cellInOrder] = currentCell;
582  cellInOrder++;
583 
584  // find if the neighbours have been visited
585  label i0 = (currentCell == 0 ? 0 : offsets[currentCell-1]);
586  label i1 = offsets[currentCell];
587 
588  for (label i = i0; i < i1; i++)
589  {
590  label nbr = cellCellAddressing.m()[i];
591 
592  if (!cellRemoved(nbr) && visited.get(nbr) == 0)
593  {
594  // not visited, add to the list
595  nextCell.append(nbr);
596  }
597  }
598  }
599  }
600  while (nextCell.size());
601  }
602  }
603 
604 
605  // Now we have new-to-old in newOrder.
606  newOrder.setSize(cellInOrder);
607 
608  // Invert to get old-to-new. Make sure removed (i.e. unmapped) cells are -1.
609  oldToNew = invert(cellCellAddressing.size(), newOrder);
610 
611  return cellInOrder;
612 }
613 
614 
615 // Determine order for faces:
616 // - upper-triangular order for internal faces
617 // - external faces after internal faces and in patch order.
618 void Foam::polyTopoChange::getFaceOrder
619 (
620  const label nActiveFaces,
621  const labelList& cellFaces,
622  const labelList& cellFaceOffsets,
623 
624  labelList& oldToNew,
625  labelList& patchSizes,
626  labelList& patchStarts
627 ) const
628 {
629  oldToNew.setSize(faceOwner_.size());
630  oldToNew = -1;
631 
632  // First unassigned face
633  label newFaceI = 0;
634 
635  forAll(cellMap_, cellI)
636  {
637  label startOfCell = cellFaceOffsets[cellI];
638  label nFaces = cellFaceOffsets[cellI+1] - startOfCell;
639 
640  // Neighbouring cells
641  SortableList<label> nbr(nFaces);
642 
643  for (label i = 0; i < nFaces; i++)
644  {
645  label faceI = cellFaces[startOfCell + i];
646 
647  label nbrCellI = faceNeighbour_[faceI];
648 
649  if (faceI >= nActiveFaces)
650  {
651  // Retired face.
652  nbr[i] = -1;
653  }
654  else if (nbrCellI != -1)
655  {
656  // Internal face. Get cell on other side.
657  if (nbrCellI == cellI)
658  {
659  nbrCellI = faceOwner_[faceI];
660  }
661 
662  if (cellI < nbrCellI)
663  {
664  // CellI is master
665  nbr[i] = nbrCellI;
666  }
667  else
668  {
669  // nbrCell is master. Let it handle this face.
670  nbr[i] = -1;
671  }
672  }
673  else
674  {
675  // External face. Do later.
676  nbr[i] = -1;
677  }
678  }
679 
680  nbr.sort();
681 
682  forAll(nbr, i)
683  {
684  if (nbr[i] != -1)
685  {
686  oldToNew[cellFaces[startOfCell + nbr.indices()[i]]] =
687  newFaceI++;
688  }
689  }
690  }
691 
692 
693  // Pick up all patch faces in patch face order.
694  patchStarts.setSize(nPatches_);
695  patchStarts = 0;
696  patchSizes.setSize(nPatches_);
697  patchSizes = 0;
698 
699  patchStarts[0] = newFaceI;
700 
701  for (label faceI = 0; faceI < nActiveFaces; faceI++)
702  {
703  if (region_[faceI] >= 0)
704  {
705  patchSizes[region_[faceI]]++;
706  }
707  }
708 
709  label faceI = patchStarts[0];
710 
711  forAll(patchStarts, patchI)
712  {
713  patchStarts[patchI] = faceI;
714  faceI += patchSizes[patchI];
715  }
716 
717  //if (debug)
718  //{
719  // Pout<< "patchSizes:" << patchSizes << nl
720  // << "patchStarts:" << patchStarts << endl;
721  //}
722 
723  labelList workPatchStarts(patchStarts);
724 
725  for (label faceI = 0; faceI < nActiveFaces; faceI++)
726  {
727  if (region_[faceI] >= 0)
728  {
729  oldToNew[faceI] = workPatchStarts[region_[faceI]]++;
730  }
731  }
732 
733  // Retired faces.
734  for (label faceI = nActiveFaces; faceI < oldToNew.size(); faceI++)
735  {
736  oldToNew[faceI] = faceI;
737  }
738 
739  // Check done all faces.
740  forAll(oldToNew, faceI)
741  {
742  if (oldToNew[faceI] == -1)
743  {
745  (
746  "polyTopoChange::getFaceOrder"
747  "(const label, const labelList&, const labelList&)"
748  " const"
749  ) << "Did not determine new position"
750  << " for face " << faceI
751  << abort(FatalError);
752  }
753  }
754 }
755 
756 
757 // Reorder and compact faces according to map.
758 void Foam::polyTopoChange::reorderCompactFaces
759 (
760  const label newSize,
761  const labelList& oldToNew
762 )
763 {
764  reorder(oldToNew, faces_);
765  faces_.setCapacity(newSize);
766 
767  reorder(oldToNew, region_);
768  region_.setCapacity(newSize);
769 
770  reorder(oldToNew, faceOwner_);
771  faceOwner_.setCapacity(newSize);
772 
773  reorder(oldToNew, faceNeighbour_);
774  faceNeighbour_.setCapacity(newSize);
775 
776  // Update faceMaps.
777  reorder(oldToNew, faceMap_);
778  faceMap_.setCapacity(newSize);
779 
780  renumberReverseMap(oldToNew, reverseFaceMap_);
781 
782  renumberKey(oldToNew, faceFromPoint_);
783  renumberKey(oldToNew, faceFromEdge_);
784  inplaceReorder(oldToNew, flipFaceFlux_);
785  flipFaceFlux_.setCapacity(newSize);
786 
787  reorder(oldToNew, faceZone_);
788  faceZone_.setCapacity(newSize);
789 
790  inplaceReorder(oldToNew, faceZoneFlip_);
791  faceZoneFlip_.setCapacity(newSize);
792 }
793 
794 
795 // Compact all and orders points and faces:
796 // - points into internal followed by external points
797 // - internalfaces upper-triangular
798 // - externalfaces after internal ones.
799 void Foam::polyTopoChange::compact
800 (
801  const bool orderCells,
802  const bool orderPoints,
803  label& nInternalPoints,
804  labelList& patchSizes,
805  labelList& patchStarts
806 )
807 {
808  points_.shrink();
809  pointMap_.shrink();
810  reversePointMap_.shrink();
811  pointZone_.shrink();
812 
813  faces_.shrink();
814  region_.shrink();
815  faceOwner_.shrink();
816  faceNeighbour_.shrink();
817  faceMap_.shrink();
818  reverseFaceMap_.shrink();
819  faceZone_.shrink();
820 
821  cellMap_.shrink();
822  reverseCellMap_.shrink();
823  cellZone_.shrink();
824 
825 
826  // Compact points
827  label nActivePoints = 0;
828  {
829  labelList localPointMap(points_.size(), -1);
830  label newPointI = 0;
831 
832  if (!orderPoints)
833  {
834  nInternalPoints = -1;
835 
836  forAll(points_, pointI)
837  {
838  if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
839  {
840  localPointMap[pointI] = newPointI++;
841  }
842  }
843  nActivePoints = newPointI;
844  }
845  else
846  {
847  forAll(points_, pointI)
848  {
849  if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
850  {
851  nActivePoints++;
852  }
853  }
854 
855  // Mark boundary points
856  forAll(faceOwner_, faceI)
857  {
858  if
859  (
860  !faceRemoved(faceI)
861  && faceOwner_[faceI] >= 0
862  && faceNeighbour_[faceI] < 0
863  )
864  {
865  // Valid boundary face
866  const face& f = faces_[faceI];
867 
868  forAll(f, fp)
869  {
870  label pointI = f[fp];
871 
872  if (localPointMap[pointI] == -1)
873  {
874  if
875  (
876  pointRemoved(pointI)
877  || retiredPoints_.found(pointI)
878  )
879  {
880  FatalErrorIn("polyTopoChange::compact(..)")
881  << "Removed or retired point " << pointI
882  << " in face " << f
883  << " at position " << faceI << endl
884  << "Probably face has not been adapted for"
885  << " removed points." << abort(FatalError);
886  }
887  localPointMap[pointI] = newPointI++;
888  }
889  }
890  }
891  }
892 
893  label nBoundaryPoints = newPointI;
894  nInternalPoints = nActivePoints - nBoundaryPoints;
895 
896  // Move the boundary addressing up
897  forAll(localPointMap, pointI)
898  {
899  if (localPointMap[pointI] != -1)
900  {
901  localPointMap[pointI] += nInternalPoints;
902  }
903  }
904 
905  newPointI = 0;
906 
907  // Mark internal points
908  forAll(faceOwner_, faceI)
909  {
910  if
911  (
912  !faceRemoved(faceI)
913  && faceOwner_[faceI] >= 0
914  && faceNeighbour_[faceI] >= 0
915  )
916  {
917  // Valid internal face
918  const face& f = faces_[faceI];
919 
920  forAll(f, fp)
921  {
922  label pointI = f[fp];
923 
924  if (localPointMap[pointI] == -1)
925  {
926  if
927  (
928  pointRemoved(pointI)
929  || retiredPoints_.found(pointI)
930  )
931  {
932  FatalErrorIn("polyTopoChange::compact(..)")
933  << "Removed or retired point " << pointI
934  << " in face " << f
935  << " at position " << faceI << endl
936  << "Probably face has not been adapted for"
937  << " removed points." << abort(FatalError);
938  }
939  localPointMap[pointI] = newPointI++;
940  }
941  }
942  }
943  }
944 
945  if (newPointI != nInternalPoints)
946  {
947  FatalErrorIn("polyTopoChange::compact(..)")
948  << "Problem." << abort(FatalError);
949  }
950  newPointI = nActivePoints;
951  }
952 
953  forAllConstIter(labelHashSet, retiredPoints_, iter)
954  {
955  localPointMap[iter.key()] = newPointI++;
956  }
957 
958 
959  if (debug)
960  {
961  Pout<< "Points : active:" << nActivePoints
962  << " removed:" << points_.size()-newPointI << endl;
963  }
964 
965  reorder(localPointMap, points_);
966  points_.setCapacity(newPointI);
967 
968  // Update pointMaps
969  reorder(localPointMap, pointMap_);
970  pointMap_.setCapacity(newPointI);
971  renumberReverseMap(localPointMap, reversePointMap_);
972 
973  reorder(localPointMap, pointZone_);
974  pointZone_.setCapacity(newPointI);
975  renumber(localPointMap, retiredPoints_);
976 
977  // Use map to relabel face vertices
978  forAll(faces_, faceI)
979  {
980  face& f = faces_[faceI];
981 
982  //labelList oldF(f);
983  renumberCompact(localPointMap, f);
984 
985  if (!faceRemoved(faceI) && f.size() < 3)
986  {
987  FatalErrorIn("polyTopoChange::compact(..)")
988  << "Created illegal face " << f
989  //<< " from face " << oldF
990  << " at position:" << faceI
991  << " when filtering removed points"
992  << abort(FatalError);
993  }
994  }
995  }
996 
997 
998  // Compact faces.
999  {
1000  labelList localFaceMap(faces_.size(), -1);
1001  label newFaceI = 0;
1002 
1003  forAll(faces_, faceI)
1004  {
1005  if (!faceRemoved(faceI) && faceOwner_[faceI] >= 0)
1006  {
1007  localFaceMap[faceI] = newFaceI++;
1008  }
1009  }
1010  nActiveFaces_ = newFaceI;
1011 
1012  forAll(faces_, faceI)
1013  {
1014  if (!faceRemoved(faceI) && faceOwner_[faceI] < 0)
1015  {
1016  // Retired face
1017  localFaceMap[faceI] = newFaceI++;
1018  }
1019  }
1020 
1021  if (debug)
1022  {
1023  Pout<< "Faces : active:" << nActiveFaces_
1024  << " removed:" << faces_.size()-newFaceI << endl;
1025  }
1026 
1027  // Reorder faces.
1028  reorderCompactFaces(newFaceI, localFaceMap);
1029  }
1030 
1031  // Compact cells.
1032  {
1033  labelList localCellMap;
1034  label newCellI;
1035 
1036  if (orderCells)
1037  {
1038  // Construct cellCell addressing
1039  CompactListList<label> cellCells;
1040  makeCellCells(nActiveFaces_, cellCells);
1041 
1042  // Cell ordering (based on bandCompression). Handles removed cells.
1043  newCellI = getCellOrder(cellCells, localCellMap);
1044  }
1045  else
1046  {
1047  // Compact out removed cells
1048  localCellMap.setSize(cellMap_.size());
1049  localCellMap = -1;
1050 
1051  newCellI = 0;
1052  forAll(cellMap_, cellI)
1053  {
1054  if (!cellRemoved(cellI))
1055  {
1056  localCellMap[cellI] = newCellI++;
1057  }
1058  }
1059  }
1060 
1061  if (debug)
1062  {
1063  Pout<< "Cells : active:" << newCellI
1064  << " removed:" << cellMap_.size()-newCellI << endl;
1065  }
1066 
1067  // Renumber -if cells reordered or -if cells removed
1068  if (orderCells || (newCellI != cellMap_.size()))
1069  {
1070  reorder(localCellMap, cellMap_);
1071  cellMap_.setCapacity(newCellI);
1072  renumberReverseMap(localCellMap, reverseCellMap_);
1073 
1074  reorder(localCellMap, cellZone_);
1075  cellZone_.setCapacity(newCellI);
1076 
1077  renumberKey(localCellMap, cellFromPoint_);
1078  renumberKey(localCellMap, cellFromEdge_);
1079  renumberKey(localCellMap, cellFromFace_);
1080 
1081  // Renumber owner/neighbour. Take into account if neighbour suddenly
1082  // gets lower cell than owner.
1083  forAll(faceOwner_, faceI)
1084  {
1085  label own = faceOwner_[faceI];
1086  label nei = faceNeighbour_[faceI];
1087 
1088  if (own >= 0)
1089  {
1090  // Update owner
1091  faceOwner_[faceI] = localCellMap[own];
1092 
1093  if (nei >= 0)
1094  {
1095  // Update neighbour.
1096  faceNeighbour_[faceI] = localCellMap[nei];
1097 
1098  // Check if face needs reversing.
1099  if
1100  (
1101  faceNeighbour_[faceI] >= 0
1102  && faceNeighbour_[faceI] < faceOwner_[faceI]
1103  )
1104  {
1105  faces_[faceI] = faces_[faceI].reverseFace();
1106  Swap(faceOwner_[faceI], faceNeighbour_[faceI]);
1107  flipFaceFlux_[faceI] =
1108  (
1109  flipFaceFlux_[faceI]
1110  ? 0
1111  : 1
1112  );
1113  faceZoneFlip_[faceI] =
1114  (
1115  faceZoneFlip_[faceI]
1116  ? 0
1117  : 1
1118  );
1119  }
1120  }
1121  }
1122  else if (nei >= 0)
1123  {
1124  // Update neighbour.
1125  faceNeighbour_[faceI] = localCellMap[nei];
1126  }
1127  }
1128  }
1129  }
1130 
1131  // Reorder faces into upper-triangular and patch ordering
1132  {
1133  // Create cells (packed storage)
1134  labelList cellFaces;
1135  labelList cellFaceOffsets;
1136  makeCells(nActiveFaces_, cellFaces, cellFaceOffsets);
1137 
1138  // Do upper triangular order and patch sorting
1139  labelList localFaceMap;
1140  getFaceOrder
1141  (
1142  nActiveFaces_,
1143  cellFaces,
1144  cellFaceOffsets,
1145 
1146  localFaceMap,
1147  patchSizes,
1148  patchStarts
1149  );
1150 
1151  // Reorder faces.
1152  reorderCompactFaces(localFaceMap.size(), localFaceMap);
1153  }
1154 }
1155 
1156 
1157 // Find faces to interpolate to create value for new face. Only used if
1158 // face was inflated from edge or point. Internal faces should only be
1159 // created from internal faces, external faces only from external faces
1160 // (and ideally the same patch)
1161 // Is bit problematic if there are no faces to select, i.e. in polyDualMesh
1162 // an internal face can be created from a boundary edge with no internal
1163 // faces connected to it.
1164 Foam::labelList Foam::polyTopoChange::selectFaces
1165 (
1166  const primitiveMesh& mesh,
1167  const labelList& faceLabels,
1168  const bool internalFacesOnly
1169 )
1170 {
1171  label nFaces = 0;
1172 
1173  forAll(faceLabels, i)
1174  {
1175  label faceI = faceLabels[i];
1176 
1177  if (internalFacesOnly == mesh.isInternalFace(faceI))
1178  {
1179  nFaces++;
1180  }
1181  }
1182 
1183  labelList collectedFaces;
1184 
1185  if (nFaces == 0)
1186  {
1187  // Did not find any faces of the correct type so just use any old
1188  // face.
1189  collectedFaces = faceLabels;
1190  }
1191  else
1192  {
1193  collectedFaces.setSize(nFaces);
1194 
1195  nFaces = 0;
1196 
1197  forAll(faceLabels, i)
1198  {
1199  label faceI = faceLabels[i];
1200 
1201  if (internalFacesOnly == mesh.isInternalFace(faceI))
1202  {
1203  collectedFaces[nFaces++] = faceI;
1204  }
1205  }
1206  }
1207 
1208  return collectedFaces;
1209 }
1210 
1211 
1212 // Calculate pointMap per patch (so from patch point label to old patch point
1213 // label)
1214 void Foam::polyTopoChange::calcPatchPointMap
1215 (
1216  const List<Map<label> >& oldPatchMeshPointMaps,
1217  const polyBoundaryMesh& boundary,
1218  labelListList& patchPointMap
1219 ) const
1220 {
1221  patchPointMap.setSize(boundary.size());
1222 
1223  forAll(boundary, patchI)
1224  {
1225  const labelList& meshPoints = boundary[patchI].meshPoints();
1226 
1227  const Map<label>& oldMeshPointMap = oldPatchMeshPointMaps[patchI];
1228 
1229  labelList& curPatchPointRnb = patchPointMap[patchI];
1230 
1231  curPatchPointRnb.setSize(meshPoints.size());
1232 
1233  forAll(meshPoints, i)
1234  {
1235  if (meshPoints[i] < pointMap_.size())
1236  {
1237  // Check if old point was part of same patch
1238  Map<label>::const_iterator ozmpmIter = oldMeshPointMap.find
1239  (
1240  pointMap_[meshPoints[i]]
1241  );
1242 
1243  if (ozmpmIter != oldMeshPointMap.end())
1244  {
1245  curPatchPointRnb[i] = ozmpmIter();
1246  }
1247  else
1248  {
1249  curPatchPointRnb[i] = -1;
1250  }
1251  }
1252  else
1253  {
1254  curPatchPointRnb[i] = -1;
1255  }
1256  }
1257  }
1258 }
1259 
1260 
1261 void Foam::polyTopoChange::calcFaceInflationMaps
1262 (
1263  const polyMesh& mesh,
1264  List<objectMap>& facesFromPoints,
1265  List<objectMap>& facesFromEdges,
1266  List<objectMap>& facesFromFaces
1267 ) const
1268 {
1269  // Faces inflated from points
1270  // ~~~~~~~~~~~~~~~~~~~~~~~~~~
1271 
1272  facesFromPoints.setSize(faceFromPoint_.size());
1273 
1274  if (faceFromPoint_.size())
1275  {
1276  label nFacesFromPoints = 0;
1277 
1278  // Collect all still existing faces connected to this point.
1279  forAllConstIter(Map<label>, faceFromPoint_, iter)
1280  {
1281  label newFaceI = iter.key();
1282 
1283  if (region_[newFaceI] == -1)
1284  {
1285  // Get internal faces using point on old mesh
1286  facesFromPoints[nFacesFromPoints++] = objectMap
1287  (
1288  newFaceI,
1289  selectFaces
1290  (
1291  mesh,
1292  mesh.pointFaces()[iter()],
1293  true
1294  )
1295  );
1296  }
1297  else
1298  {
1299  // Get patch faces using point on old mesh
1300  facesFromPoints[nFacesFromPoints++] = objectMap
1301  (
1302  newFaceI,
1303  selectFaces
1304  (
1305  mesh,
1306  mesh.pointFaces()[iter()],
1307  false
1308  )
1309  );
1310  }
1311  }
1312  }
1313 
1314 
1315  // Faces inflated from edges
1316  // ~~~~~~~~~~~~~~~~~~~~~~~~~
1317 
1318  facesFromEdges.setSize(faceFromEdge_.size());
1319 
1320  if (faceFromEdge_.size())
1321  {
1322  label nFacesFromEdges = 0;
1323 
1324  // Collect all still existing faces connected to this edge.
1325  forAllConstIter(Map<label>, faceFromEdge_, iter)
1326  {
1327  label newFaceI = iter.key();
1328 
1329  if (region_[newFaceI] == -1)
1330  {
1331  // Get internal faces using edge on old mesh
1332  facesFromEdges[nFacesFromEdges++] = objectMap
1333  (
1334  newFaceI,
1335  selectFaces
1336  (
1337  mesh,
1338  mesh.edgeFaces(iter()),
1339  true
1340  )
1341  );
1342  }
1343  else
1344  {
1345  // Get patch faces using edge on old mesh
1346  facesFromEdges[nFacesFromEdges++] = objectMap
1347  (
1348  newFaceI,
1349  selectFaces
1350  (
1351  mesh,
1352  mesh.edgeFaces(iter()),
1353  false
1354  )
1355  );
1356  }
1357  }
1358  }
1359 
1360 
1361  // Faces from face merging
1362  // ~~~~~~~~~~~~~~~~~~~~~~~
1363 
1364  getMergeSets
1365  (
1366  reverseFaceMap_,
1367  faceMap_,
1368  facesFromFaces
1369  );
1370 }
1371 
1372 
1373 void Foam::polyTopoChange::calcCellInflationMaps
1374 (
1375  const polyMesh& mesh,
1376  List<objectMap>& cellsFromPoints,
1377  List<objectMap>& cellsFromEdges,
1378  List<objectMap>& cellsFromFaces,
1379  List<objectMap>& cellsFromCells
1380 ) const
1381 {
1382  cellsFromPoints.setSize(cellFromPoint_.size());
1383 
1384  if (cellFromPoint_.size())
1385  {
1386  label nCellsFromPoints = 0;
1387 
1388  // Collect all still existing faces connected to this point.
1389  forAllConstIter(Map<label>, cellFromPoint_, iter)
1390  {
1391  cellsFromPoints[nCellsFromPoints++] = objectMap
1392  (
1393  iter.key(),
1394  mesh.pointCells()[iter()]
1395  );
1396  }
1397  }
1398 
1399 
1400  cellsFromEdges.setSize(cellFromEdge_.size());
1401 
1402  if (cellFromEdge_.size())
1403  {
1404  label nCellsFromEdges = 0;
1405 
1406  // Collect all still existing faces connected to this point.
1407  forAllConstIter(Map<label>, cellFromEdge_, iter)
1408  {
1409  cellsFromEdges[nCellsFromEdges++] = objectMap
1410  (
1411  iter.key(),
1412  mesh.edgeCells()[iter()]
1413  );
1414  }
1415  }
1416 
1417 
1418  cellsFromFaces.setSize(cellFromFace_.size());
1419 
1420  if (cellFromFace_.size())
1421  {
1422  label nCellsFromFaces = 0;
1423 
1424  labelList twoCells(2);
1425 
1426  // Collect all still existing faces connected to this point.
1427  forAllConstIter(Map<label>, cellFromFace_, iter)
1428  {
1429  label oldFaceI = iter();
1430 
1431  if (mesh.isInternalFace(oldFaceI))
1432  {
1433  twoCells[0] = mesh.faceOwner()[oldFaceI];
1434  twoCells[1] = mesh.faceNeighbour()[oldFaceI];
1435  cellsFromFaces[nCellsFromFaces++] = objectMap
1436  (
1437  iter.key(),
1438  twoCells
1439  );
1440  }
1441  else
1442  {
1443  cellsFromFaces[nCellsFromFaces++] = objectMap
1444  (
1445  iter.key(),
1446  labelList(1, mesh.faceOwner()[oldFaceI])
1447  );
1448  }
1449  }
1450  }
1451 
1452 
1453  // Cells from cell merging
1454  // ~~~~~~~~~~~~~~~~~~~~~~~
1455 
1456  getMergeSets
1457  (
1458  reverseCellMap_,
1459  cellMap_,
1460  cellsFromCells
1461  );
1462 }
1463 
1464 
1465 void Foam::polyTopoChange::resetZones
1466 (
1467  const polyMesh& mesh,
1468  polyMesh& newMesh,
1469  labelListList& pointZoneMap,
1470  labelListList& faceZoneFaceMap,
1471  labelListList& cellZoneMap
1472 ) const
1473 {
1474  // pointZones
1475  // ~~~~~~~~~~
1476 
1477  pointZoneMap.setSize(mesh.pointZones().size());
1478  {
1479  const pointZoneMesh& pointZones = mesh.pointZones();
1480 
1481  // Count points per zone
1482 
1483  labelList nPoints(pointZones.size(), 0);
1484 
1485  forAll(pointZone_, pointI)
1486  {
1487  label zoneI = pointZone_[pointI];
1488 
1489  if (zoneI >= pointZones.size())
1490  {
1491  FatalErrorIn
1492  (
1493  "resetZones(const polyMesh&, polyMesh&, labelListList&"
1494  "labelListList&, labelListList&)"
1495  ) << "Illegal zoneID " << zoneI << " for point "
1496  << pointI << " coord " << mesh.points()[pointI]
1497  << abort(FatalError);
1498  }
1499 
1500  if (zoneI >= 0)
1501  {
1502  nPoints[zoneI]++;
1503  }
1504  }
1505 
1506  // Distribute points per zone
1507 
1508  labelListList addressing(pointZones.size());
1509  forAll(addressing, zoneI)
1510  {
1511  addressing[zoneI].setSize(nPoints[zoneI]);
1512  }
1513  nPoints = 0;
1514 
1515  forAll(pointZone_, pointI)
1516  {
1517  label zoneI = pointZone_[pointI];
1518  if (zoneI >= 0)
1519  {
1520  addressing[zoneI][nPoints[zoneI]++] = pointI;
1521  }
1522  }
1523  // Sort the addressing
1524  forAll(addressing, zoneI)
1525  {
1526  stableSort(addressing[zoneI]);
1527  }
1528 
1529  // So now we both have old zones and the new addressing.
1530  // Invert the addressing to get pointZoneMap.
1531  forAll(addressing, zoneI)
1532  {
1533  const pointZone& oldZone = pointZones[zoneI];
1534  const labelList& newZoneAddr = addressing[zoneI];
1535 
1536  labelList& curPzRnb = pointZoneMap[zoneI];
1537  curPzRnb.setSize(newZoneAddr.size());
1538 
1539  forAll(newZoneAddr, i)
1540  {
1541  if (newZoneAddr[i] < pointMap_.size())
1542  {
1543  curPzRnb[i] = oldZone.whichPoint(pointMap_[newZoneAddr[i]]);
1544  }
1545  else
1546  {
1547  curPzRnb[i] = -1;
1548  }
1549  }
1550  }
1551 
1552  // Reset the addresing on the zone
1553  newMesh.pointZones().clearAddressing();
1554  forAll(newMesh.pointZones(), zoneI)
1555  {
1556  if (debug)
1557  {
1558  Pout<< "pointZone:" << zoneI
1559  << " name:" << newMesh.pointZones()[zoneI].name()
1560  << " size:" << addressing[zoneI].size()
1561  << endl;
1562  }
1563 
1564  newMesh.pointZones()[zoneI] = addressing[zoneI];
1565  }
1566  }
1567 
1568 
1569  // faceZones
1570  // ~~~~~~~~~
1571 
1572  faceZoneFaceMap.setSize(mesh.faceZones().size());
1573  {
1574  const faceZoneMesh& faceZones = mesh.faceZones();
1575 
1576  labelList nFaces(faceZones.size(), 0);
1577 
1578  forAll(faceZone_, faceI)
1579  {
1580  label zoneI = faceZone_[faceI];
1581 
1582  if (zoneI >= faceZones.size())
1583  {
1584  FatalErrorIn
1585  (
1586  "resetZones(const polyMesh&, polyMesh&, labelListList&"
1587  "labelListList&, labelListList&)"
1588  ) << "Illegal zoneID " << zoneI << " for face "
1589  << faceI
1590  << abort(FatalError);
1591  }
1592  if (zoneI >= 0)
1593  {
1594  nFaces[zoneI]++;
1595  }
1596  }
1597 
1598  labelListList addressing(faceZones.size());
1599  boolListList flipMode(faceZones.size());
1600 
1601  forAll(addressing, zoneI)
1602  {
1603  addressing[zoneI].setSize(nFaces[zoneI]);
1604  flipMode[zoneI].setSize(nFaces[zoneI]);
1605  }
1606  nFaces = 0;
1607 
1608  forAll(faceZone_, faceI)
1609  {
1610  label zoneI = faceZone_[faceI];
1611 
1612  if (zoneI >= 0)
1613  {
1614  label index = nFaces[zoneI]++;
1615  addressing[zoneI][index] = faceI;
1616  flipMode[zoneI][index] = faceZoneFlip_[faceI];
1617  }
1618  }
1619  // Sort the addressing
1620  forAll(addressing, zoneI)
1621  {
1622  labelList newToOld;
1623  sortedOrder(addressing[zoneI], newToOld);
1624  {
1625  labelList newAddressing(addressing[zoneI].size());
1626  forAll(newAddressing, i)
1627  {
1628  newAddressing[i] = addressing[zoneI][newToOld[i]];
1629  }
1630  addressing[zoneI].transfer(newAddressing);
1631  }
1632  {
1633  boolList newFlipMode(flipMode[zoneI].size());
1634  forAll(newFlipMode, i)
1635  {
1636  newFlipMode[i] = flipMode[zoneI][newToOld[i]];
1637  }
1638  flipMode[zoneI].transfer(newFlipMode);
1639  }
1640  }
1641 
1642  // So now we both have old zones and the new addressing.
1643  // Invert the addressing to get faceZoneFaceMap.
1644  forAll(addressing, zoneI)
1645  {
1646  const faceZone& oldZone = faceZones[zoneI];
1647  const labelList& newZoneAddr = addressing[zoneI];
1648 
1649  labelList& curFzFaceRnb = faceZoneFaceMap[zoneI];
1650 
1651  curFzFaceRnb.setSize(newZoneAddr.size());
1652 
1653  forAll(newZoneAddr, i)
1654  {
1655  if (newZoneAddr[i] < faceMap_.size())
1656  {
1657  curFzFaceRnb[i] =
1658  oldZone.whichFace(faceMap_[newZoneAddr[i]]);
1659  }
1660  else
1661  {
1662  curFzFaceRnb[i] = -1;
1663  }
1664  }
1665  }
1666 
1667 
1668  // Reset the addresing on the zone
1669  newMesh.faceZones().clearAddressing();
1670  forAll(newMesh.faceZones(), zoneI)
1671  {
1672  if (debug)
1673  {
1674  Pout<< "faceZone:" << zoneI
1675  << " name:" << newMesh.faceZones()[zoneI].name()
1676  << " size:" << addressing[zoneI].size()
1677  << endl;
1678  }
1679 
1680  newMesh.faceZones()[zoneI].resetAddressing
1681  (
1682  addressing[zoneI],
1683  flipMode[zoneI]
1684  );
1685  }
1686  }
1687 
1688 
1689  // cellZones
1690  // ~~~~~~~~~
1691 
1692  cellZoneMap.setSize(mesh.cellZones().size());
1693  {
1694  const cellZoneMesh& cellZones = mesh.cellZones();
1695 
1696  labelList nCells(cellZones.size(), 0);
1697 
1698  forAll(cellZone_, cellI)
1699  {
1700  label zoneI = cellZone_[cellI];
1701 
1702  if (zoneI >= cellZones.size())
1703  {
1704  FatalErrorIn
1705  (
1706  "resetZones(const polyMesh&, polyMesh&, labelListList&"
1707  "labelListList&, labelListList&)"
1708  ) << "Illegal zoneID " << zoneI << " for cell "
1709  << cellI << abort(FatalError);
1710  }
1711 
1712  if (zoneI >= 0)
1713  {
1714  nCells[zoneI]++;
1715  }
1716  }
1717 
1718  labelListList addressing(cellZones.size());
1719  forAll(addressing, zoneI)
1720  {
1721  addressing[zoneI].setSize(nCells[zoneI]);
1722  }
1723  nCells = 0;
1724 
1725  forAll(cellZone_, cellI)
1726  {
1727  label zoneI = cellZone_[cellI];
1728 
1729  if (zoneI >= 0)
1730  {
1731  addressing[zoneI][nCells[zoneI]++] = cellI;
1732  }
1733  }
1734  // Sort the addressing
1735  forAll(addressing, zoneI)
1736  {
1737  stableSort(addressing[zoneI]);
1738  }
1739 
1740  // So now we both have old zones and the new addressing.
1741  // Invert the addressing to get cellZoneMap.
1742  forAll(addressing, zoneI)
1743  {
1744  const cellZone& oldZone = cellZones[zoneI];
1745  const labelList& newZoneAddr = addressing[zoneI];
1746 
1747  labelList& curCellRnb = cellZoneMap[zoneI];
1748 
1749  curCellRnb.setSize(newZoneAddr.size());
1750 
1751  forAll(newZoneAddr, i)
1752  {
1753  if (newZoneAddr[i] < cellMap_.size())
1754  {
1755  curCellRnb[i] =
1756  oldZone.whichCell(cellMap_[newZoneAddr[i]]);
1757  }
1758  else
1759  {
1760  curCellRnb[i] = -1;
1761  }
1762  }
1763  }
1764 
1765  // Reset the addresing on the zone
1766  newMesh.cellZones().clearAddressing();
1767  forAll(newMesh.cellZones(), zoneI)
1768  {
1769  if (debug)
1770  {
1771  Pout<< "cellZone:" << zoneI
1772  << " name:" << newMesh.cellZones()[zoneI].name()
1773  << " size:" << addressing[zoneI].size()
1774  << endl;
1775  }
1776 
1777  newMesh.cellZones()[zoneI] = addressing[zoneI];
1778  }
1779  }
1780 }
1781 
1782 
1783 void Foam::polyTopoChange::calcFaceZonePointMap
1784 (
1785  const polyMesh& mesh,
1786  const List<Map<label> >& oldFaceZoneMeshPointMaps,
1787  labelListList& faceZonePointMap
1788 ) const
1789 {
1790  const faceZoneMesh& faceZones = mesh.faceZones();
1791 
1792  faceZonePointMap.setSize(faceZones.size());
1793 
1794  forAll(faceZones, zoneI)
1795  {
1796  const faceZone& newZone = faceZones[zoneI];
1797 
1798  const labelList& newZoneMeshPoints = newZone().meshPoints();
1799 
1800  const Map<label>& oldZoneMeshPointMap = oldFaceZoneMeshPointMaps[zoneI];
1801 
1802  labelList& curFzPointRnb = faceZonePointMap[zoneI];
1803 
1804  curFzPointRnb.setSize(newZoneMeshPoints.size());
1805 
1806  forAll(newZoneMeshPoints, pointI)
1807  {
1808  if (newZoneMeshPoints[pointI] < pointMap_.size())
1809  {
1810  Map<label>::const_iterator ozmpmIter =
1811  oldZoneMeshPointMap.find
1812  (
1813  pointMap_[newZoneMeshPoints[pointI]]
1814  );
1815 
1816  if (ozmpmIter != oldZoneMeshPointMap.end())
1817  {
1818  curFzPointRnb[pointI] = ozmpmIter();
1819  }
1820  else
1821  {
1822  curFzPointRnb[pointI] = -1;
1823  }
1824  }
1825  else
1826  {
1827  curFzPointRnb[pointI] = -1;
1828  }
1829  }
1830  }
1831 }
1832 
1833 
1834 Foam::face Foam::polyTopoChange::rotateFace
1835 (
1836  const face& f,
1837  const label nPos
1838 )
1839 {
1840  face newF(f.size());
1841 
1842  forAll(f, fp)
1843  {
1844  label fp1 = (fp + nPos) % f.size();
1845 
1846  if (fp1 < 0)
1847  {
1848  fp1 += f.size();
1849  }
1850 
1851  newF[fp1] = f[fp];
1852  }
1853 
1854  return newF;
1855 }
1856 
1857 
1858 void Foam::polyTopoChange::reorderCoupledFaces
1859 (
1860  const bool syncParallel,
1861  const polyBoundaryMesh& boundary,
1862  const labelList& patchStarts,
1863  const labelList& patchSizes,
1864  const pointField& points
1865 )
1866 {
1867  // Mapping for faces (old to new). Extends over all mesh faces for
1868  // convenience (could be just the external faces)
1869  labelList oldToNew(identity(faces_.size()));
1870 
1871  // Rotation on new faces.
1872  labelList rotation(faces_.size(), 0);
1873 
1874  // Send ordering
1875  forAll(boundary, patchI)
1876  {
1877  if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1878  {
1879  boundary[patchI].initOrder
1880  (
1882  (
1883  SubList<face>
1884  (
1885  faces_,
1886  patchSizes[patchI],
1887  patchStarts[patchI]
1888  ),
1889  points
1890  )
1891  );
1892  }
1893  }
1894 
1895  // Receive and calculate ordering
1896 
1897  bool anyChanged = false;
1898 
1899  forAll(boundary, patchI)
1900  {
1901  if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1902  {
1903  labelList patchFaceMap(patchSizes[patchI], -1);
1904  labelList patchFaceRotation(patchSizes[patchI], 0);
1905 
1906  bool changed = boundary[patchI].order
1907  (
1909  (
1910  SubList<face>
1911  (
1912  faces_,
1913  patchSizes[patchI],
1914  patchStarts[patchI]
1915  ),
1916  points
1917  ),
1918  patchFaceMap,
1919  patchFaceRotation
1920  );
1921 
1922  if (changed)
1923  {
1924  // Merge patch face reordering into mesh face reordering table
1925  label start = patchStarts[patchI];
1926 
1927  forAll(patchFaceMap, patchFaceI)
1928  {
1929  oldToNew[patchFaceI + start] =
1930  start + patchFaceMap[patchFaceI];
1931  }
1932 
1933  forAll(patchFaceRotation, patchFaceI)
1934  {
1935  rotation[patchFaceI + start] =
1936  patchFaceRotation[patchFaceI];
1937  }
1938 
1939  anyChanged = true;
1940  }
1941  }
1942  }
1943 
1944  if (syncParallel)
1945  {
1946  reduce(anyChanged, orOp<bool>());
1947  }
1948 
1949  if (anyChanged)
1950  {
1951  // Reorder faces according to oldToNew.
1952  reorderCompactFaces(oldToNew.size(), oldToNew);
1953 
1954  // Rotate faces (rotation is already in new face indices).
1955  forAll(rotation, faceI)
1956  {
1957  if (rotation[faceI] != 0)
1958  {
1959  faces_[faceI] = rotateFace(faces_[faceI], rotation[faceI]);
1960  }
1961  }
1962  }
1963 }
1964 
1965 
1966 void Foam::polyTopoChange::compactAndReorder
1967 (
1968  const polyMesh& mesh,
1969  const bool syncParallel,
1970  const bool orderCells,
1971  const bool orderPoints,
1972 
1973  label& nInternalPoints,
1974  pointField& newPoints,
1975  labelList& patchSizes,
1976  labelList& patchStarts,
1977  List<objectMap>& pointsFromPoints,
1978  List<objectMap>& facesFromPoints,
1979  List<objectMap>& facesFromEdges,
1980  List<objectMap>& facesFromFaces,
1981  List<objectMap>& cellsFromPoints,
1982  List<objectMap>& cellsFromEdges,
1983  List<objectMap>& cellsFromFaces,
1984  List<objectMap>& cellsFromCells,
1985  List<Map<label> >& oldPatchMeshPointMaps,
1986  labelList& oldPatchNMeshPoints,
1987  labelList& oldPatchStarts,
1988  List<Map<label> >& oldFaceZoneMeshPointMaps
1989 )
1990 {
1991  if (mesh.boundaryMesh().size() != nPatches_)
1992  {
1993  FatalErrorIn("polyTopoChange::compactAndReorder(..)")
1994  << "polyTopoChange was constructed with a mesh with "
1995  << nPatches_ << " patches." << endl
1996  << "The mesh now provided has a different number of patches "
1997  << mesh.boundaryMesh().size()
1998  << " which is illegal" << endl
1999  << abort(FatalError);
2000  }
2001 
2002  // Remove any holes from points/faces/cells and sort faces.
2003  // Sets nActiveFaces_.
2004  compact(orderCells, orderPoints, nInternalPoints, patchSizes, patchStarts);
2005 
2006  // Transfer points to pointField. points_ are now cleared!
2007  // Only done since e.g. reorderCoupledFaces requires pointField.
2008  newPoints.transfer(points_);
2009 
2010  // Reorder any coupled faces
2011  reorderCoupledFaces
2012  (
2013  syncParallel,
2014  mesh.boundaryMesh(),
2015  patchStarts,
2016  patchSizes,
2017  newPoints
2018  );
2019 
2020 
2021  // Calculate inflation/merging maps
2022  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2023  // These are for the new face(/point/cell) the old faces whose value
2024  // needs to be
2025  // averaged/summed to get the new value. These old faces come either from
2026  // merged old faces (face remove into other face),
2027  // the old edgeFaces (inflate from edge) or the old pointFaces (inflate
2028  // from point). As an additional complexity will use only internal faces
2029  // to create new value for internal face and vice versa only patch
2030  // faces to to create patch face value.
2031 
2032  // For point only point merging
2033  getMergeSets
2034  (
2035  reversePointMap_,
2036  pointMap_,
2037  pointsFromPoints
2038  );
2039 
2040  calcFaceInflationMaps
2041  (
2042  mesh,
2043  facesFromPoints,
2044  facesFromEdges,
2045  facesFromFaces
2046  );
2047 
2048  calcCellInflationMaps
2049  (
2050  mesh,
2051  cellsFromPoints,
2052  cellsFromEdges,
2053  cellsFromFaces,
2054  cellsFromCells
2055  );
2056 
2057  // Clear inflation info
2058  {
2059  faceFromPoint_.clearStorage();
2060  faceFromEdge_.clearStorage();
2061 
2062  cellFromPoint_.clearStorage();
2063  cellFromEdge_.clearStorage();
2064  cellFromFace_.clearStorage();
2065  }
2066 
2067 
2068  const polyBoundaryMesh& boundary = mesh.boundaryMesh();
2069 
2070  // Grab patch mesh point maps
2071  oldPatchMeshPointMaps.setSize(boundary.size());
2072  oldPatchNMeshPoints.setSize(boundary.size());
2073  oldPatchStarts.setSize(boundary.size());
2074 
2075  forAll(boundary, patchI)
2076  {
2077  // Copy old face zone mesh point maps
2078  oldPatchMeshPointMaps[patchI] = boundary[patchI].meshPointMap();
2079  oldPatchNMeshPoints[patchI] = boundary[patchI].meshPoints().size();
2080  oldPatchStarts[patchI] = boundary[patchI].start();
2081  }
2082 
2083  // Grab old face zone mesh point maps.
2084  // These need to be saved before resetting the mesh and are used
2085  // later on to calculate the faceZone pointMaps.
2086  oldFaceZoneMeshPointMaps.setSize(mesh.faceZones().size());
2087 
2088  forAll(mesh.faceZones(), zoneI)
2089  {
2090  const faceZone& oldZone = mesh.faceZones()[zoneI];
2091 
2092  oldFaceZoneMeshPointMaps[zoneI] = oldZone().meshPointMap();
2093  }
2094 }
2095 
2096 
2097 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
2098 
2099 // Construct from components
2100 Foam::polyTopoChange::polyTopoChange(const label nPatches, const bool strict)
2101 :
2102  strict_(strict),
2103  nPatches_(nPatches),
2104  points_(0),
2105  pointMap_(0),
2106  reversePointMap_(0),
2107  pointZone_(0),
2108  retiredPoints_(0),
2109  faces_(0),
2110  region_(0),
2111  faceOwner_(0),
2112  faceNeighbour_(0),
2113  faceMap_(0),
2114  reverseFaceMap_(0),
2115  faceFromPoint_(0),
2116  faceFromEdge_(0),
2117  flipFaceFlux_(0),
2118  faceZone_(0),
2119  faceZoneFlip_(0),
2120  nActiveFaces_(0),
2121  cellMap_(0),
2122  reverseCellMap_(0),
2123  cellFromPoint_(0),
2124  cellFromEdge_(0),
2125  cellFromFace_(0),
2126  cellZone_(0)
2127 {}
2128 
2129 
2130 // Construct from components
2133  const polyMesh& mesh,
2134  const bool strict
2135 )
2136 :
2137  strict_(strict),
2138  nPatches_(0),
2139  points_(0),
2140  pointMap_(0),
2141  reversePointMap_(0),
2142  pointZone_(0),
2143  retiredPoints_(0),
2144  faces_(0),
2145  region_(0),
2146  faceOwner_(0),
2147  faceNeighbour_(0),
2148  faceMap_(0),
2149  reverseFaceMap_(0),
2150  faceFromPoint_(0),
2151  faceFromEdge_(0),
2152  flipFaceFlux_(0),
2153  faceZone_(0),
2154  faceZoneFlip_(0),
2155  nActiveFaces_(0),
2156  cellMap_(0),
2157  reverseCellMap_(0),
2158  cellFromPoint_(0),
2159  cellFromEdge_(0),
2160  cellFromFace_(0),
2161  cellZone_(0)
2162 {
2163  addMesh
2164  (
2165  mesh,
2166  identity(mesh.boundaryMesh().size()),
2167  identity(mesh.pointZones().size()),
2168  identity(mesh.faceZones().size()),
2169  identity(mesh.cellZones().size())
2170  );
2171 }
2172 
2173 
2174 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2175 
2177 {
2178  points_.clearStorage();
2179  pointMap_.clearStorage();
2180  reversePointMap_.clearStorage();
2181  pointZone_.clearStorage();
2182  retiredPoints_.clearStorage();
2183 
2184  faces_.clearStorage();
2185  region_.clearStorage();
2186  faceOwner_.clearStorage();
2187  faceNeighbour_.clearStorage();
2188  faceMap_.clearStorage();
2189  reverseFaceMap_.clearStorage();
2190  faceFromPoint_.clearStorage();
2191  faceFromEdge_.clearStorage();
2192  flipFaceFlux_.clearStorage();
2193  faceZone_.clearStorage();
2194  faceZoneFlip_.clearStorage();
2195  nActiveFaces_ = 0;
2196 
2197  cellMap_.clearStorage();
2198  reverseCellMap_.clearStorage();
2199  cellZone_.clearStorage();
2200  cellFromPoint_.clearStorage();
2201  cellFromEdge_.clearStorage();
2202  cellFromFace_.clearStorage();
2203 }
2204 
2205 
2208  const polyMesh& mesh,
2209  const labelList& patchMap,
2210  const labelList& pointZoneMap,
2211  const labelList& faceZoneMap,
2212  const labelList& cellZoneMap
2213 )
2214 {
2215  label maxRegion = nPatches_ - 1;
2216  forAll(patchMap, i)
2217  {
2218  maxRegion = max(maxRegion, patchMap[i]);
2219  }
2220  nPatches_ = maxRegion + 1;
2221 
2222 
2223  // Add points
2224  {
2225  const pointField& points = mesh.points();
2226  const pointZoneMesh& pointZones = mesh.pointZones();
2227 
2228  // Extend
2229  points_.setCapacity(points_.size() + points.size());
2230  pointMap_.setCapacity(pointMap_.size() + points.size());
2231  reversePointMap_.setCapacity(reversePointMap_.size() + points.size());
2232  pointZone_.setCapacity(pointZone_.size() + points.size());
2233 
2234  // Precalc offset zones
2235  labelList newZoneID(points.size(), -1);
2236 
2237  forAll(pointZones, zoneI)
2238  {
2239  const labelList& pointLabels = pointZones[zoneI];
2240 
2241  forAll(pointLabels, j)
2242  {
2243  newZoneID[pointLabels[j]] = pointZoneMap[zoneI];
2244  }
2245  }
2246 
2247  // Add points in mesh order
2248  for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
2249  {
2250  addPoint
2251  (
2252  points[pointI],
2253  pointI,
2254  newZoneID[pointI],
2255  true
2256  );
2257  }
2258  }
2259 
2260  // Add cells
2261  {
2262  const cellZoneMesh& cellZones = mesh.cellZones();
2263 
2264  // Resize
2265 
2266  // Note: polyMesh does not allow retired cells anymore. So allCells
2267  // always equals nCells
2268  label nAllCells = mesh.nCells();
2269 
2270  cellMap_.setCapacity(cellMap_.size() + nAllCells);
2271  reverseCellMap_.setCapacity(reverseCellMap_.size() + nAllCells);
2272  cellFromPoint_.resize(cellFromPoint_.size() + nAllCells/100);
2273  cellFromEdge_.resize(cellFromEdge_.size() + nAllCells/100);
2274  cellFromFace_.resize(cellFromFace_.size() + nAllCells/100);
2275  cellZone_.setCapacity(cellZone_.size() + nAllCells);
2276 
2277 
2278  // Precalc offset zones
2279  labelList newZoneID(nAllCells, -1);
2280 
2281  forAll(cellZones, zoneI)
2282  {
2283  const labelList& cellLabels = cellZones[zoneI];
2284 
2285  forAll(cellLabels, j)
2286  {
2287  label cellI = cellLabels[j];
2288 
2289  if (newZoneID[cellI] != -1)
2290  {
2291  WarningIn
2292  (
2293  "polyTopoChange::addMesh"
2294  "(const polyMesh&, const labelList&,"
2295  "const labelList&, const labelList&,"
2296  "const labelList&)"
2297  ) << "Cell:" << cellI
2298  << " centre:" << mesh.cellCentres()[cellI]
2299  << " is in two zones:"
2300  << cellZones[newZoneID[cellI]].name()
2301  << " and " << cellZones[zoneI].name() << endl
2302  << " This is not supported."
2303  << " Continuing with first zone only." << endl;
2304  }
2305  else
2306  {
2307  newZoneID[cellI] = cellZoneMap[zoneI];
2308  }
2309  }
2310  }
2311 
2312  // Add cells in mesh order
2313  for (label cellI = 0; cellI < nAllCells; cellI++)
2314  {
2315  // Add cell from cell
2316  addCell(-1, -1, -1, cellI, newZoneID[cellI]);
2317  }
2318  }
2319 
2320  // Add faces
2321  {
2322  const polyBoundaryMesh& patches = mesh.boundaryMesh();
2323  const faceList& faces = mesh.faces();
2324  const labelList& faceOwner = mesh.faceOwner();
2325  const labelList& faceNeighbour = mesh.faceNeighbour();
2326  const faceZoneMesh& faceZones = mesh.faceZones();
2327 
2328  // Resize
2329  label nAllFaces = mesh.faces().size();
2330 
2331  faces_.setCapacity(faces_.size() + nAllFaces);
2332  region_.setCapacity(region_.size() + nAllFaces);
2333  faceOwner_.setCapacity(faceOwner_.size() + nAllFaces);
2334  faceNeighbour_.setCapacity(faceNeighbour_.size() + nAllFaces);
2335  faceMap_.setCapacity(faceMap_.size() + nAllFaces);
2336  reverseFaceMap_.setCapacity(reverseFaceMap_.size() + nAllFaces);
2337  faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/100);
2338  faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/100);
2339  flipFaceFlux_.setCapacity(flipFaceFlux_.size() + nAllFaces);
2340  faceZone_.setCapacity(faceZone_.size() + nAllFaces);
2341  faceZoneFlip_.setCapacity(faceZoneFlip_.size() + nAllFaces);
2342 
2343 
2344  // Precalc offset zones
2345  labelList newZoneID(nAllFaces, -1);
2346  boolList zoneFlip(nAllFaces, false);
2347 
2348  forAll(faceZones, zoneI)
2349  {
2350  const labelList& faceLabels = faceZones[zoneI];
2351  const boolList& flipMap = faceZones[zoneI].flipMap();
2352 
2353  forAll(faceLabels, j)
2354  {
2355  newZoneID[faceLabels[j]] = faceZoneMap[zoneI];
2356  zoneFlip[faceLabels[j]] = flipMap[j];
2357  }
2358  }
2359 
2360  // Add faces in mesh order
2361 
2362  // 1. Internal faces
2363  for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
2364  {
2365  addFace
2366  (
2367  faces[faceI],
2368  faceOwner[faceI],
2369  faceNeighbour[faceI],
2370  -1, // masterPointID
2371  -1, // masterEdgeID
2372  faceI, // masterFaceID
2373  false, // flipFaceFlux
2374  -1, // patchID
2375  newZoneID[faceI], // zoneID
2376  zoneFlip[faceI] // zoneFlip
2377  );
2378  }
2379 
2380  // 2. Patch faces
2381  forAll(patches, patchI)
2382  {
2383  const polyPatch& pp = patches[patchI];
2384 
2385  if (pp.start() != faces_.size())
2386  {
2387  FatalErrorIn
2388  (
2389  "polyTopoChange::polyTopoChange"
2390  "(const polyMesh& mesh, const bool strict)"
2391  ) << "Problem : "
2392  << "Patch " << pp.name() << " starts at " << pp.start()
2393  << endl
2394  << "Current face counter at " << faces_.size() << endl
2395  << "Are patches in incremental order?"
2396  << abort(FatalError);
2397  }
2398  forAll(pp, patchFaceI)
2399  {
2400  label faceI = pp.start() + patchFaceI;
2401 
2402  addFace
2403  (
2404  faces[faceI],
2405  faceOwner[faceI],
2406  -1, // neighbour
2407  -1, // masterPointID
2408  -1, // masterEdgeID
2409  faceI, // masterFaceID
2410  false, // flipFaceFlux
2411  patchMap[patchI], // patchID
2412  newZoneID[faceI], // zoneID
2413  zoneFlip[faceI] // zoneFlip
2414  );
2415  }
2416  }
2417  }
2418 }
2419 
2420 
2423  const label nPoints,
2424  const label nFaces,
2425  const label nCells
2426 )
2427 {
2428  points_.setCapacity(nPoints);
2429  pointMap_.setCapacity(nPoints);
2430  reversePointMap_.setCapacity(nPoints);
2431  pointZone_.setCapacity(nPoints);
2432 
2433  faces_.setCapacity(nFaces);
2434  region_.setCapacity(nFaces);
2435  faceOwner_.setCapacity(nFaces);
2436  faceNeighbour_.setCapacity(nFaces);
2437  faceMap_.setCapacity(nFaces);
2438  reverseFaceMap_.setCapacity(nFaces);
2439  faceFromPoint_.resize(faceFromPoint_.size() + nFaces/100);
2440  faceFromEdge_.resize(faceFromEdge_.size() + nFaces/100);
2441  flipFaceFlux_.setCapacity(nFaces);
2442  faceZone_.setCapacity(nFaces);
2443  faceZoneFlip_.setCapacity(nFaces);
2444 
2445  cellMap_.setCapacity(nCells);
2446  reverseCellMap_.setCapacity(nCells);
2447  cellFromPoint_.resize(cellFromPoint_.size() + nCells/100);
2448  cellFromEdge_.resize(cellFromEdge_.size() + nCells/100);
2449  cellFromFace_.resize(cellFromFace_.size() + nCells/100);
2450  cellZone_.setCapacity(nCells);
2451 }
2452 
2453 
2455 {
2456  if (isType<polyAddPoint>(action))
2457  {
2458  const polyAddPoint& pap = refCast<const polyAddPoint>(action);
2459 
2460  return addPoint
2461  (
2462  pap.newPoint(),
2463  pap.masterPointID(),
2464  pap.zoneID(),
2465  pap.inCell()
2466  );
2467  }
2468  else if (isType<polyModifyPoint>(action))
2469  {
2470  const polyModifyPoint& pmp = refCast<const polyModifyPoint>(action);
2471 
2472  modifyPoint
2473  (
2474  pmp.pointID(),
2475  pmp.newPoint(),
2476  pmp.zoneID(),
2477  pmp.inCell()
2478  );
2479 
2480  return -1;
2481  }
2482  else if (isType<polyRemovePoint>(action))
2483  {
2484  const polyRemovePoint& prp = refCast<const polyRemovePoint>(action);
2485 
2486  removePoint(prp.pointID(), prp.mergePointID());
2487 
2488  return -1;
2489  }
2490  else if (isType<polyAddFace>(action))
2491  {
2492  const polyAddFace& paf = refCast<const polyAddFace>(action);
2493 
2494  return addFace
2495  (
2496  paf.newFace(),
2497  paf.owner(),
2498  paf.neighbour(),
2499  paf.masterPointID(),
2500  paf.masterEdgeID(),
2501  paf.masterFaceID(),
2502  paf.flipFaceFlux(),
2503  paf.patchID(),
2504  paf.zoneID(),
2505  paf.zoneFlip()
2506  );
2507  }
2508  else if (isType<polyModifyFace>(action))
2509  {
2510  const polyModifyFace& pmf = refCast<const polyModifyFace>(action);
2511 
2512  modifyFace
2513  (
2514  pmf.newFace(),
2515  pmf.faceID(),
2516  pmf.owner(),
2517  pmf.neighbour(),
2518  pmf.flipFaceFlux(),
2519  pmf.patchID(),
2520  pmf.zoneID(),
2521  pmf.zoneFlip()
2522  );
2523 
2524  return -1;
2525  }
2526  else if (isType<polyRemoveFace>(action))
2527  {
2528  const polyRemoveFace& prf = refCast<const polyRemoveFace>(action);
2529 
2530  removeFace(prf.faceID(), prf.mergeFaceID());
2531 
2532  return -1;
2533  }
2534  else if (isType<polyAddCell>(action))
2535  {
2536  const polyAddCell& pac = refCast<const polyAddCell>(action);
2537 
2538  return addCell
2539  (
2540  pac.masterPointID(),
2541  pac.masterEdgeID(),
2542  pac.masterFaceID(),
2543  pac.masterCellID(),
2544  pac.zoneID()
2545  );
2546  }
2547  else if (isType<polyModifyCell>(action))
2548  {
2549  const polyModifyCell& pmc = refCast<const polyModifyCell>(action);
2550 
2551  if (pmc.removeFromZone())
2552  {
2553  modifyCell(pmc.cellID(), -1);
2554  }
2555  else
2556  {
2557  modifyCell(pmc.cellID(), pmc.zoneID());
2558  }
2559 
2560  return -1;
2561  }
2562  else if (isType<polyRemoveCell>(action))
2563  {
2564  const polyRemoveCell& prc = refCast<const polyRemoveCell>(action);
2565 
2566  removeCell(prc.cellID(), prc.mergeCellID());
2567 
2568  return -1;
2569  }
2570  else
2571  {
2572  FatalErrorIn
2573  (
2574  "label polyTopoChange::setAction(const topoAction& action)"
2575  ) << "Unknown type of topoChange: " << action.type()
2576  << abort(FatalError);
2577 
2578  // Dummy return to keep compiler happy
2579  return -1;
2580  }
2581 }
2582 
2583 
2586  const point& pt,
2587  const label masterPointID,
2588  const label zoneID,
2589  const bool inCell
2590 )
2591 {
2592  label pointI = points_.size();
2593 
2594  points_.append(pt);
2595  pointMap_.append(masterPointID);
2596  reversePointMap_.append(pointI);
2597  pointZone_.append(zoneID);
2598 
2599  if (!inCell)
2600  {
2601  retiredPoints_.insert(pointI);
2602  }
2603 
2604  return pointI;
2605 }
2606 
2607 
2610  const label pointI,
2611  const point& pt,
2612  const label newZoneID,
2613  const bool inCell
2614 )
2615 {
2616  if (pointI < 0 || pointI >= points_.size())
2617  {
2618  FatalErrorIn
2619  (
2620  "polyTopoChange::modifyPoint(const label, const point&)"
2621  ) << "illegal point label " << pointI << endl
2622  << "Valid point labels are 0 .. " << points_.size()-1
2623  << abort(FatalError);
2624  }
2625  if (pointRemoved(pointI) || pointMap_[pointI] == -1)
2626  {
2627  FatalErrorIn
2628  (
2629  "polyTopoChange::modifyPoint(const label, const point&)"
2630  ) << "point " << pointI << " already marked for removal"
2631  << abort(FatalError);
2632  }
2633  points_[pointI] = pt;
2634  pointZone_[pointI] = newZoneID;
2635 
2636  if (inCell)
2637  {
2638  retiredPoints_.erase(pointI);
2639  }
2640  else
2641  {
2642  retiredPoints_.insert(pointI);
2643  }
2644 }
2645 
2646 
2648 {
2649  if (newPoints.size() != points_.size())
2650  {
2651  FatalErrorIn("polyTopoChange::movePoints(const pointField&)")
2652  << "illegal pointField size." << endl
2653  << "Size:" << newPoints.size() << endl
2654  << "Points in mesh:" << points_.size()
2655  << abort(FatalError);
2656  }
2657 
2658  forAll(points_, pointI)
2659  {
2660  points_[pointI] = newPoints[pointI];
2661  }
2662 }
2663 
2664 
2667  const label pointI,
2668  const label mergePointI
2669 )
2670 {
2671  if (pointI < 0 || pointI >= points_.size())
2672  {
2673  FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2674  << "illegal point label " << pointI << endl
2675  << "Valid point labels are 0 .. " << points_.size()-1
2676  << abort(FatalError);
2677  }
2678 
2679  if
2680  (
2681  strict_
2682  && (pointRemoved(pointI) || pointMap_[pointI] == -1)
2683  )
2684  {
2685  FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2686  << "point " << pointI << " already marked for removal" << nl
2687  << "Point:" << points_[pointI] << " pointMap:" << pointMap_[pointI]
2688  << abort(FatalError);
2689  }
2690 
2691  if (pointI == mergePointI)
2692  {
2693  FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2694  << "Cannot remove/merge point " << pointI << " onto itself."
2695  << abort(FatalError);
2696  }
2697 
2698  points_[pointI] = greatPoint;
2699  pointMap_[pointI] = -1;
2700  if (mergePointI >= 0)
2701  {
2702  reversePointMap_[pointI] = -mergePointI-2;
2703  }
2704  else
2705  {
2706  reversePointMap_[pointI] = -1;
2707  }
2708  pointZone_[pointI] = -1;
2709  retiredPoints_.erase(pointI);
2710 }
2711 
2712 
2715  const face& f,
2716  const label own,
2717  const label nei,
2718  const label masterPointID,
2719  const label masterEdgeID,
2720  const label masterFaceID,
2721  const bool flipFaceFlux,
2722  const label patchID,
2723  const label zoneID,
2724  const bool zoneFlip
2725 )
2726 {
2727  // Check validity
2728  if (debug)
2729  {
2730  checkFace(f, -1, own, nei, patchID, zoneID);
2731  }
2732 
2733  label faceI = faces_.size();
2734 
2735  faces_.append(f);
2736  region_.append(patchID);
2737  faceOwner_.append(own);
2738  faceNeighbour_.append(nei);
2739 
2740  if (masterPointID >= 0)
2741  {
2742  faceMap_.append(-1);
2743  faceFromPoint_.insert(faceI, masterPointID);
2744  }
2745  else if (masterEdgeID >= 0)
2746  {
2747  faceMap_.append(-1);
2748  faceFromEdge_.insert(faceI, masterEdgeID);
2749  }
2750  else if (masterFaceID >= 0)
2751  {
2752  faceMap_.append(masterFaceID);
2753  }
2754  else
2755  {
2756  // Allow inflate-from-nothing?
2757  //FatalErrorIn("polyTopoChange::addFace")
2758  // << "Need to specify a master point, edge or face"
2759  // << "face:" << f << " own:" << own << " nei:" << nei
2760  // << abort(FatalError);
2761  faceMap_.append(-1);
2762  }
2763  reverseFaceMap_.append(faceI);
2764 
2765  flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2766  faceZone_.append(zoneID);
2767  faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2768 
2769  return faceI;
2770 }
2771 
2772 
2775  const face& f,
2776  const label faceI,
2777  const label own,
2778  const label nei,
2779  const bool flipFaceFlux,
2780  const label patchID,
2781  const label zoneID,
2782  const bool zoneFlip
2783 )
2784 {
2785  // Check validity
2786  if (debug)
2787  {
2788  checkFace(f, faceI, own, nei, patchID, zoneID);
2789  }
2790 
2791  faces_[faceI] = f;
2792  faceOwner_[faceI] = own;
2793  faceNeighbour_[faceI] = nei;
2794  region_[faceI] = patchID;
2795 
2796  flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2797  faceZone_[faceI] = zoneID;
2798  faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2799 }
2800 
2801 
2802 void Foam::polyTopoChange::removeFace(const label faceI, const label mergeFaceI)
2803 {
2804  if (faceI < 0 || faceI >= faces_.size())
2805  {
2806  FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2807  << "illegal face label " << faceI << endl
2808  << "Valid face labels are 0 .. " << faces_.size()-1
2809  << abort(FatalError);
2810  }
2811 
2812  if
2813  (
2814  strict_
2815  && (faceRemoved(faceI) || faceMap_[faceI] == -1)
2816  )
2817  {
2818  FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2819  << "face " << faceI
2820  << " already marked for removal"
2821  << abort(FatalError);
2822  }
2823 
2824  faces_[faceI].setSize(0);
2825  region_[faceI] = -1;
2826  faceOwner_[faceI] = -1;
2827  faceNeighbour_[faceI] = -1;
2828  faceMap_[faceI] = -1;
2829  if (mergeFaceI >= 0)
2830  {
2831  reverseFaceMap_[faceI] = -mergeFaceI-2;
2832  }
2833  else
2834  {
2835  reverseFaceMap_[faceI] = -1;
2836  }
2837  faceFromEdge_.erase(faceI);
2838  faceFromPoint_.erase(faceI);
2839  flipFaceFlux_[faceI] = 0;
2840  faceZone_[faceI] = -1;
2841  faceZoneFlip_[faceI] = 0;
2842 }
2843 
2844 
2847  const label masterPointID,
2848  const label masterEdgeID,
2849  const label masterFaceID,
2850  const label masterCellID,
2851  const label zoneID
2852 )
2853 {
2854  label cellI = cellMap_.size();
2855 
2856  if (masterPointID >= 0)
2857  {
2858  cellMap_.append(-1);
2859  cellFromPoint_.insert(cellI, masterPointID);
2860  }
2861  else if (masterEdgeID >= 0)
2862  {
2863  cellMap_.append(-1);
2864  cellFromEdge_.insert(cellI, masterEdgeID);
2865  }
2866  else if (masterFaceID >= 0)
2867  {
2868  cellMap_.append(-1);
2869  cellFromFace_.insert(cellI, masterFaceID);
2870  }
2871  else
2872  {
2873  cellMap_.append(masterCellID);
2874  }
2875  reverseCellMap_.append(cellI);
2876  cellZone_.append(zoneID);
2877 
2878  return cellI;
2879 }
2880 
2881 
2884  const label cellI,
2885  const label zoneID
2886 )
2887 {
2888  cellZone_[cellI] = zoneID;
2889 }
2890 
2891 
2892 void Foam::polyTopoChange::removeCell(const label cellI, const label mergeCellI)
2893 {
2894  if (cellI < 0 || cellI >= cellMap_.size())
2895  {
2896  FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2897  << "illegal cell label " << cellI << endl
2898  << "Valid cell labels are 0 .. " << cellMap_.size()-1
2899  << abort(FatalError);
2900  }
2901 
2902  if (strict_ && cellMap_[cellI] == -2)
2903  {
2904  FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2905  << "cell " << cellI
2906  << " already marked for removal"
2907  << abort(FatalError);
2908  }
2909 
2910  cellMap_[cellI] = -2;
2911  if (mergeCellI >= 0)
2912  {
2913  reverseCellMap_[cellI] = -mergeCellI-2;
2914  }
2915  else
2916  {
2917  reverseCellMap_[cellI] = -1;
2918  }
2919  cellFromPoint_.erase(cellI);
2920  cellFromEdge_.erase(cellI);
2921  cellFromFace_.erase(cellI);
2922  cellZone_[cellI] = -1;
2923 }
2924 
2925 
2928  polyMesh& mesh,
2929  const bool inflate,
2930  const bool syncParallel,
2931  const bool orderCells,
2932  const bool orderPoints
2933 )
2934 {
2935  if (debug)
2936  {
2937  Pout<< "polyTopoChange::changeMesh"
2938  << "(polyMesh&, const bool, const bool, const bool, const bool)"
2939  << endl;
2940  }
2941 
2942  if (debug)
2943  {
2944  Pout<< "Old mesh:" << nl;
2945  writeMeshStats(mesh, Pout);
2946  }
2947 
2948  // new mesh points
2949  pointField newPoints;
2950  // number of internal points
2951  label nInternalPoints;
2952  // patch slicing
2953  labelList patchSizes;
2954  labelList patchStarts;
2955  // inflate maps
2956  List<objectMap> pointsFromPoints;
2957  List<objectMap> facesFromPoints;
2958  List<objectMap> facesFromEdges;
2959  List<objectMap> facesFromFaces;
2960  List<objectMap> cellsFromPoints;
2961  List<objectMap> cellsFromEdges;
2962  List<objectMap> cellsFromFaces;
2963  List<objectMap> cellsFromCells;
2964  // old mesh info
2965  List<Map<label> > oldPatchMeshPointMaps;
2966  labelList oldPatchNMeshPoints;
2967  labelList oldPatchStarts;
2968  List<Map<label> > oldFaceZoneMeshPointMaps;
2969 
2970  // Compact, reorder patch faces and calculate mesh/patch maps.
2971  compactAndReorder
2972  (
2973  mesh,
2974  syncParallel,
2975  orderCells,
2976  orderPoints,
2977 
2978  nInternalPoints,
2979  newPoints,
2980  patchSizes,
2981  patchStarts,
2982  pointsFromPoints,
2983  facesFromPoints,
2984  facesFromEdges,
2985  facesFromFaces,
2986  cellsFromPoints,
2987  cellsFromEdges,
2988  cellsFromFaces,
2989  cellsFromCells,
2990  oldPatchMeshPointMaps,
2991  oldPatchNMeshPoints,
2992  oldPatchStarts,
2993  oldFaceZoneMeshPointMaps
2994  );
2995 
2996  const label nOldPoints(mesh.nPoints());
2997  const label nOldFaces(mesh.nFaces());
2998  const label nOldCells(mesh.nCells());
2999 
3000 
3001  // Change the mesh
3002  // ~~~~~~~~~~~~~~~
3003  // This will invalidate any addressing so better make sure you have
3004  // all the information you need!!!
3005 
3006  if (inflate)
3007  {
3008  // Keep (renumbered) mesh points, store new points in map for inflation
3009  // (appended points (i.e. from nowhere) get value zero)
3010  pointField renumberedMeshPoints(newPoints.size());
3011 
3012  forAll(pointMap_, newPointI)
3013  {
3014  label oldPointI = pointMap_[newPointI];
3015 
3016  if (oldPointI >= 0)
3017  {
3018  renumberedMeshPoints[newPointI] = mesh.points()[oldPointI];
3019  }
3020  else
3021  {
3022  renumberedMeshPoints[newPointI] = vector::zero;
3023  }
3024  }
3025 
3026  mesh.resetPrimitives
3027  (
3028  xferMove(renumberedMeshPoints),
3029  faces_.xfer(),
3030  faceOwner_.xfer(),
3031  faceNeighbour_.xfer(),
3032  patchSizes,
3033  patchStarts,
3034  syncParallel
3035  );
3036 
3037  mesh.changing(true);
3038  }
3039  else
3040  {
3041  // Set new points.
3042  mesh.resetPrimitives
3043  (
3044  xferMove(newPoints),
3045  faces_.xfer(),
3046  faceOwner_.xfer(),
3047  faceNeighbour_.xfer(),
3048  patchSizes,
3049  patchStarts,
3050  syncParallel
3051  );
3052  mesh.changing(true);
3053  }
3054 
3055  // Clear out primitives
3056  {
3057  retiredPoints_.clearStorage();
3058  region_.clearStorage();
3059  }
3060 
3061 
3062  if (debug)
3063  {
3064  // Some stats on changes
3065  label nAdd, nInflate, nMerge, nRemove;
3066  countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3067  Pout<< "Points:"
3068  << " added(from point):" << nAdd
3069  << " added(from nothing):" << nInflate
3070  << " merged(into other point):" << nMerge
3071  << " removed:" << nRemove
3072  << nl;
3073 
3074  countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3075  Pout<< "Faces:"
3076  << " added(from face):" << nAdd
3077  << " added(inflated):" << nInflate
3078  << " merged(into other face):" << nMerge
3079  << " removed:" << nRemove
3080  << nl;
3081 
3082  countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3083  Pout<< "Cells:"
3084  << " added(from cell):" << nAdd
3085  << " added(inflated):" << nInflate
3086  << " merged(into other cell):" << nMerge
3087  << " removed:" << nRemove
3088  << nl
3089  << endl;
3090  }
3091 
3092  if (debug)
3093  {
3094  Pout<< "New mesh:" << nl;
3095  writeMeshStats(mesh, Pout);
3096  }
3097 
3098 
3099  // Zones
3100  // ~~~~~
3101 
3102  // Inverse of point/face/cell zone addressing.
3103  // For every preserved point/face/cells in zone give the old position.
3104  // For added points, the index is set to -1
3105  labelListList pointZoneMap(mesh.pointZones().size());
3106  labelListList faceZoneFaceMap(mesh.faceZones().size());
3107  labelListList cellZoneMap(mesh.cellZones().size());
3108 
3109  resetZones(mesh, mesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3110 
3111  // Clear zone info
3112  {
3113  pointZone_.clearStorage();
3114  faceZone_.clearStorage();
3115  faceZoneFlip_.clearStorage();
3116  cellZone_.clearStorage();
3117  }
3118 
3119 
3120  // Patch point renumbering
3121  // For every preserved point on a patch give the old position.
3122  // For added points, the index is set to -1
3123  labelListList patchPointMap(mesh.boundaryMesh().size());
3124  calcPatchPointMap
3125  (
3126  oldPatchMeshPointMaps,
3127  mesh.boundaryMesh(),
3128  patchPointMap
3129  );
3130 
3131  // Create the face zone mesh point renumbering
3132  labelListList faceZonePointMap(mesh.faceZones().size());
3133  calcFaceZonePointMap(mesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3134 
3135  labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3136 
3137  return autoPtr<mapPolyMesh>
3138  (
3139  new mapPolyMesh
3140  (
3141  mesh,
3142  nOldPoints,
3143  nOldFaces,
3144  nOldCells,
3145 
3146  pointMap_,
3147  pointsFromPoints,
3148 
3149  faceMap_,
3150  facesFromPoints,
3151  facesFromEdges,
3152  facesFromFaces,
3153 
3154  cellMap_,
3155  cellsFromPoints,
3156  cellsFromEdges,
3157  cellsFromFaces,
3158  cellsFromCells,
3159 
3160  reversePointMap_,
3161  reverseFaceMap_,
3162  reverseCellMap_,
3163 
3164  flipFaceFluxSet,
3165 
3166  patchPointMap,
3167 
3168  pointZoneMap,
3169 
3170  faceZonePointMap,
3171  faceZoneFaceMap,
3172  cellZoneMap,
3173 
3174  newPoints, // if empty signals no inflation.
3175  oldPatchStarts,
3176  oldPatchNMeshPoints,
3177  true // steal storage.
3178  )
3179  );
3180 
3181  // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3182  // be invalid.
3183 }
3184 
3185 
3188  autoPtr<fvMesh>& newMeshPtr,
3189  const IOobject& io,
3190  const fvMesh& mesh,
3191  const bool syncParallel,
3192  const bool orderCells,
3193  const bool orderPoints
3194 )
3195 {
3196  if (debug)
3197  {
3198  Pout<< "polyTopoChange::changeMesh"
3199  << "(autoPtr<fvMesh>&, const IOobject&, const fvMesh&"
3200  << ", const bool, const bool, const bool)"
3201  << endl;
3202  }
3203 
3204  if (debug)
3205  {
3206  Pout<< "Old mesh:" << nl;
3207  writeMeshStats(mesh, Pout);
3208  }
3209 
3210  // new mesh points
3211  pointField newPoints;
3212  // number of internal points
3213  label nInternalPoints;
3214  // patch slicing
3215  labelList patchSizes;
3216  labelList patchStarts;
3217  // inflate maps
3218  List<objectMap> pointsFromPoints;
3219  List<objectMap> facesFromPoints;
3220  List<objectMap> facesFromEdges;
3221  List<objectMap> facesFromFaces;
3222  List<objectMap> cellsFromPoints;
3223  List<objectMap> cellsFromEdges;
3224  List<objectMap> cellsFromFaces;
3225  List<objectMap> cellsFromCells;
3226 
3227  // old mesh info
3228  List<Map<label> > oldPatchMeshPointMaps;
3229  labelList oldPatchNMeshPoints;
3230  labelList oldPatchStarts;
3231  List<Map<label> > oldFaceZoneMeshPointMaps;
3232 
3233  // Compact, reorder patch faces and calculate mesh/patch maps.
3234  compactAndReorder
3235  (
3236  mesh,
3237  syncParallel,
3238  orderCells,
3239  orderPoints,
3240 
3241  nInternalPoints,
3242  newPoints,
3243  patchSizes,
3244  patchStarts,
3245  pointsFromPoints,
3246  facesFromPoints,
3247  facesFromEdges,
3248  facesFromFaces,
3249  cellsFromPoints,
3250  cellsFromEdges,
3251  cellsFromFaces,
3252  cellsFromCells,
3253  oldPatchMeshPointMaps,
3254  oldPatchNMeshPoints,
3255  oldPatchStarts,
3256  oldFaceZoneMeshPointMaps
3257  );
3258 
3259  const label nOldPoints(mesh.nPoints());
3260  const label nOldFaces(mesh.nFaces());
3261  const label nOldCells(mesh.nCells());
3262 
3263 
3264  // Create the mesh
3265  // ~~~~~~~~~~~~~~~
3266 
3267  newMeshPtr.reset
3268  (
3269  new fvMesh
3270  (
3271  io,
3272  xferMove(newPoints),
3273  faces_.xfer(),
3274  faceOwner_.xfer(),
3275  faceNeighbour_.xfer()
3276  )
3277  );
3278  fvMesh& newMesh = newMeshPtr();
3279 
3280  // Clear out primitives
3281  {
3282  retiredPoints_.clearStorage();
3283  region_.clearStorage();
3284  }
3285 
3286 
3287  if (debug)
3288  {
3289  // Some stats on changes
3290  label nAdd, nInflate, nMerge, nRemove;
3291  countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3292  Pout<< "Points:"
3293  << " added(from point):" << nAdd
3294  << " added(from nothing):" << nInflate
3295  << " merged(into other point):" << nMerge
3296  << " removed:" << nRemove
3297  << nl;
3298 
3299  countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3300  Pout<< "Faces:"
3301  << " added(from face):" << nAdd
3302  << " added(inflated):" << nInflate
3303  << " merged(into other face):" << nMerge
3304  << " removed:" << nRemove
3305  << nl;
3306 
3307  countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3308  Pout<< "Cells:"
3309  << " added(from cell):" << nAdd
3310  << " added(inflated):" << nInflate
3311  << " merged(into other cell):" << nMerge
3312  << " removed:" << nRemove
3313  << nl
3314  << endl;
3315  }
3316 
3317 
3318  {
3319  const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
3320 
3321  List<polyPatch*> newBoundary(oldPatches.size());
3322 
3323  forAll(oldPatches, patchI)
3324  {
3325  newBoundary[patchI] = oldPatches[patchI].clone
3326  (
3327  newMesh.boundaryMesh(),
3328  patchI,
3329  patchSizes[patchI],
3330  patchStarts[patchI]
3331  ).ptr();
3332  }
3333  newMesh.addFvPatches(newBoundary);
3334  }
3335 
3336 
3337  // Zones
3338  // ~~~~~
3339 
3340  // Start off from empty zones.
3341  const pointZoneMesh& oldPointZones = mesh.pointZones();
3342  List<pointZone*> pZonePtrs(oldPointZones.size());
3343  {
3344  forAll(oldPointZones, i)
3345  {
3346  pZonePtrs[i] = new pointZone
3347  (
3348  oldPointZones[i].name(),
3349  labelList(0),
3350  i,
3351  newMesh.pointZones()
3352  );
3353  }
3354  }
3355 
3356  const faceZoneMesh& oldFaceZones = mesh.faceZones();
3357  List<faceZone*> fZonePtrs(oldFaceZones.size());
3358  {
3359  forAll(oldFaceZones, i)
3360  {
3361  fZonePtrs[i] = new faceZone
3362  (
3363  oldFaceZones[i].name(),
3364  labelList(0),
3365  boolList(0),
3366  i,
3367  newMesh.faceZones()
3368  );
3369  }
3370  }
3371 
3372  const cellZoneMesh& oldCellZones = mesh.cellZones();
3373  List<cellZone*> cZonePtrs(oldCellZones.size());
3374  {
3375  forAll(oldCellZones, i)
3376  {
3377  cZonePtrs[i] = new cellZone
3378  (
3379  oldCellZones[i].name(),
3380  labelList(0),
3381  i,
3382  newMesh.cellZones()
3383  );
3384  }
3385  }
3386 
3387  newMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
3388 
3389  // Inverse of point/face/cell zone addressing.
3390  // For every preserved point/face/cells in zone give the old position.
3391  // For added points, the index is set to -1
3392  labelListList pointZoneMap(mesh.pointZones().size());
3393  labelListList faceZoneFaceMap(mesh.faceZones().size());
3394  labelListList cellZoneMap(mesh.cellZones().size());
3395 
3396  resetZones(mesh, newMesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3397 
3398  // Clear zone info
3399  {
3400  pointZone_.clearStorage();
3401  faceZone_.clearStorage();
3402  faceZoneFlip_.clearStorage();
3403  cellZone_.clearStorage();
3404  }
3405 
3406  // Patch point renumbering
3407  // For every preserved point on a patch give the old position.
3408  // For added points, the index is set to -1
3409  labelListList patchPointMap(newMesh.boundaryMesh().size());
3410  calcPatchPointMap
3411  (
3412  oldPatchMeshPointMaps,
3413  newMesh.boundaryMesh(),
3414  patchPointMap
3415  );
3416 
3417  // Create the face zone mesh point renumbering
3418  labelListList faceZonePointMap(newMesh.faceZones().size());
3419  calcFaceZonePointMap(newMesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3420 
3421  if (debug)
3422  {
3423  Pout<< "New mesh:" << nl;
3424  writeMeshStats(mesh, Pout);
3425  }
3426 
3427  labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3428 
3429  return autoPtr<mapPolyMesh>
3430  (
3431  new mapPolyMesh
3432  (
3433  newMesh,
3434  nOldPoints,
3435  nOldFaces,
3436  nOldCells,
3437 
3438  pointMap_,
3439  pointsFromPoints,
3440 
3441  faceMap_,
3442  facesFromPoints,
3443  facesFromEdges,
3444  facesFromFaces,
3445 
3446  cellMap_,
3447  cellsFromPoints,
3448  cellsFromEdges,
3449  cellsFromFaces,
3450  cellsFromCells,
3451 
3452  reversePointMap_,
3453  reverseFaceMap_,
3454  reverseCellMap_,
3455 
3456  flipFaceFluxSet,
3457 
3458  patchPointMap,
3459 
3460  pointZoneMap,
3461 
3462  faceZonePointMap,
3463  faceZoneFaceMap,
3464  cellZoneMap,
3465 
3466  newPoints, // if empty signals no inflation.
3467  oldPatchStarts,
3468  oldPatchNMeshPoints,
3469  true // steal storage.
3470  )
3471  );
3472 
3473  // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3474  // be invalid.
3475 }
3476 
3477 
3478 // ************************ vim: set sw=4 sts=4 et: ************************ //