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