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