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-2019 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 minimize 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 pointZoneMesh& 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 faceZoneMesh& 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 cellZoneMesh& 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>>& oldFaceZoneMeshPointMaps,
1872  labelListList& faceZonePointMap
1873 ) const
1874 {
1875  const faceZoneMesh& 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& newZoneMeshPoints = newZone().meshPoints();
1884 
1885  const Map<label>& oldZoneMeshPointMap = oldFaceZoneMeshPointMaps[zoneI];
1886 
1887  labelList& curFzPointRnb = faceZonePointMap[zoneI];
1888 
1889  curFzPointRnb.setSize(newZoneMeshPoints.size());
1890 
1891  forAll(newZoneMeshPoints, pointi)
1892  {
1893  if (newZoneMeshPoints[pointi] < pointMap_.size())
1894  {
1895  Map<label>::const_iterator ozmpmIter =
1896  oldZoneMeshPointMap.find
1897  (
1898  pointMap_[newZoneMeshPoints[pointi]]
1899  );
1900 
1901  if (ozmpmIter != oldZoneMeshPointMap.end())
1902  {
1903  curFzPointRnb[pointi] = ozmpmIter();
1904  }
1905  else
1906  {
1907  curFzPointRnb[pointi] = -1;
1908  }
1909  }
1910  else
1911  {
1912  curFzPointRnb[pointi] = -1;
1913  }
1914  }
1915  }
1916 }
1917 
1918 
1919 void Foam::polyTopoChange::reorderCoupledFaces
1920 (
1921  const bool syncParallel,
1922  const polyBoundaryMesh& boundary,
1923  const labelList& patchStarts,
1924  const labelList& patchSizes,
1925  const pointField& points
1926 )
1927 {
1928  // Mapping for faces (old to new). Extends over all mesh faces for
1929  // convenience (could be just the external faces)
1930  labelList oldToNew(identity(faces_.size()));
1931 
1932  // Rotation on new faces.
1933  labelList rotation(faces_.size(), 0);
1934 
1935  PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
1936 
1937  // Send ordering
1938  forAll(boundary, patchi)
1939  {
1940  if (syncParallel || !isA<processorPolyPatch>(boundary[patchi]))
1941  {
1942  boundary[patchi].initOrder
1943  (
1944  pBufs,
1946  (
1947  SubList<face>
1948  (
1949  faces_,
1950  patchSizes[patchi],
1951  patchStarts[patchi]
1952  ),
1953  points
1954  )
1955  );
1956  }
1957  }
1958 
1959  if (syncParallel)
1960  {
1961  pBufs.finishedSends();
1962  }
1963 
1964  // Receive and calculate ordering
1965 
1966  bool anyChanged = false;
1967 
1968  forAll(boundary, patchi)
1969  {
1970  if (syncParallel || !isA<processorPolyPatch>(boundary[patchi]))
1971  {
1972  labelList patchFaceMap(patchSizes[patchi], -1);
1973  labelList patchFaceRotation(patchSizes[patchi], 0);
1974 
1975  bool changed = boundary[patchi].order
1976  (
1977  pBufs,
1979  (
1980  SubList<face>
1981  (
1982  faces_,
1983  patchSizes[patchi],
1984  patchStarts[patchi]
1985  ),
1986  points
1987  ),
1988  patchFaceMap,
1989  patchFaceRotation
1990  );
1991 
1992  if (changed)
1993  {
1994  // Merge patch face reordering into mesh face reordering table
1995  label start = patchStarts[patchi];
1996 
1997  forAll(patchFaceMap, patchFacei)
1998  {
1999  oldToNew[patchFacei + start] =
2000  start + patchFaceMap[patchFacei];
2001  }
2002 
2003  forAll(patchFaceRotation, patchFacei)
2004  {
2005  rotation[patchFacei + start] =
2006  patchFaceRotation[patchFacei];
2007  }
2008 
2009  anyChanged = true;
2010  }
2011  }
2012  }
2013 
2014  if (syncParallel)
2015  {
2016  reduce(anyChanged, orOp<bool>());
2017  }
2018 
2019  if (anyChanged)
2020  {
2021  // Reorder faces according to oldToNew.
2022  reorderCompactFaces(oldToNew.size(), oldToNew);
2023 
2024  // Rotate faces (rotation is already in new face indices).
2025  forAll(rotation, facei)
2026  {
2027  if (rotation[facei] != 0)
2028  {
2029  inplaceRotateList<List, label>(faces_[facei], rotation[facei]);
2030  }
2031  }
2032  }
2033 }
2034 
2035 
2036 void Foam::polyTopoChange::compactAndReorder
2037 (
2038  const polyMesh& mesh,
2039  const bool syncParallel,
2040  const bool orderCells,
2041  const bool orderPoints,
2042 
2043  label& nInternalPoints,
2044  pointField& newPoints,
2045  labelList& patchSizes,
2046  labelList& patchStarts,
2047  List<objectMap>& pointsFromPoints,
2048  List<objectMap>& facesFromPoints,
2049  List<objectMap>& facesFromEdges,
2050  List<objectMap>& facesFromFaces,
2051  List<objectMap>& cellsFromPoints,
2052  List<objectMap>& cellsFromEdges,
2053  List<objectMap>& cellsFromFaces,
2054  List<objectMap>& cellsFromCells,
2055  List<Map<label>>& oldPatchMeshPointMaps,
2056  labelList& oldPatchNMeshPoints,
2057  labelList& oldPatchStarts,
2058  List<Map<label>>& oldFaceZoneMeshPointMaps
2059 )
2060 {
2061  if (mesh.boundaryMesh().size() != nPatches_)
2062  {
2064  << "polyTopoChange was constructed with a mesh with "
2065  << nPatches_ << " patches." << endl
2066  << "The mesh now provided has a different number of patches "
2067  << mesh.boundaryMesh().size()
2068  << " which is illegal" << endl
2069  << abort(FatalError);
2070  }
2071 
2072  // Remove any holes from points/faces/cells and sort faces.
2073  // Sets nActiveFaces_.
2074  compact(orderCells, orderPoints, nInternalPoints, patchSizes, patchStarts);
2075 
2076  // Transfer points to pointField. points_ are now cleared!
2077  // Only done since e.g. reorderCoupledFaces requires pointField.
2078  newPoints.transfer(points_);
2079 
2080  // Reorder any coupled faces
2081  reorderCoupledFaces
2082  (
2083  syncParallel,
2084  mesh.boundaryMesh(),
2085  patchStarts,
2086  patchSizes,
2087  newPoints
2088  );
2089 
2090 
2091  // Calculate inflation/merging maps
2092  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2093  // These are for the new face(/point/cell) the old faces whose value
2094  // needs to be
2095  // averaged/summed to get the new value. These old faces come either from
2096  // merged old faces (face remove into other face),
2097  // the old edgeFaces (inflate from edge) or the old pointFaces (inflate
2098  // from point). As an additional complexity will use only internal faces
2099  // to create new value for internal face and vice versa only patch
2100  // faces to to create patch face value.
2101 
2102  // For point only point merging
2103  getMergeSets
2104  (
2105  reversePointMap_,
2106  pointMap_,
2107  pointsFromPoints
2108  );
2109 
2110  calcFaceInflationMaps
2111  (
2112  mesh,
2113  facesFromPoints,
2114  facesFromEdges,
2115  facesFromFaces
2116  );
2117 
2118  calcCellInflationMaps
2119  (
2120  mesh,
2121  cellsFromPoints,
2122  cellsFromEdges,
2123  cellsFromFaces,
2124  cellsFromCells
2125  );
2126 
2127  // Clear inflation info
2128  {
2129  faceFromPoint_.clearStorage();
2130  faceFromEdge_.clearStorage();
2131 
2132  cellFromPoint_.clearStorage();
2133  cellFromEdge_.clearStorage();
2134  cellFromFace_.clearStorage();
2135  }
2136 
2137 
2138  const polyBoundaryMesh& boundary = mesh.boundaryMesh();
2139 
2140  // Grab patch mesh point maps
2141  oldPatchMeshPointMaps.setSize(boundary.size());
2142  oldPatchNMeshPoints.setSize(boundary.size());
2143  oldPatchStarts.setSize(boundary.size());
2144 
2145  forAll(boundary, patchi)
2146  {
2147  // Copy old face zone mesh point maps
2148  oldPatchMeshPointMaps[patchi] = boundary[patchi].meshPointMap();
2149  oldPatchNMeshPoints[patchi] = boundary[patchi].meshPoints().size();
2150  oldPatchStarts[patchi] = boundary[patchi].start();
2151  }
2152 
2153  // Grab old face zone mesh point maps.
2154  // These need to be saved before resetting the mesh and are used
2155  // later on to calculate the faceZone pointMaps.
2156  oldFaceZoneMeshPointMaps.setSize(mesh.faceZones().size());
2157 
2158  forAll(mesh.faceZones(), zoneI)
2159  {
2160  const faceZone& oldZone = mesh.faceZones()[zoneI];
2161 
2162  oldFaceZoneMeshPointMaps[zoneI] = oldZone().meshPointMap();
2163  }
2164 }
2165 
2166 
2167 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
2168 
2169 // Construct from components
2171 :
2172  strict_(strict),
2173  nPatches_(nPatches),
2174  points_(0),
2175  pointMap_(0),
2176  reversePointMap_(0),
2177  pointZone_(0),
2178  retiredPoints_(0),
2179  oldPoints_(0),
2180  faces_(0),
2181  region_(0),
2182  faceOwner_(0),
2183  faceNeighbour_(0),
2184  faceMap_(0),
2185  reverseFaceMap_(0),
2186  faceFromPoint_(0),
2187  faceFromEdge_(0),
2188  flipFaceFlux_(0),
2189  faceZone_(0),
2190  faceZoneFlip_(0),
2191  nActiveFaces_(0),
2192  cellMap_(0),
2193  reverseCellMap_(0),
2194  cellFromPoint_(0),
2195  cellFromEdge_(0),
2196  cellFromFace_(0),
2197  cellZone_(0)
2198 {}
2199 
2200 
2201 // Construct from components
2204  const polyMesh& mesh,
2205  const bool strict
2206 )
2207 :
2208  strict_(strict),
2209  nPatches_(0),
2210  points_(0),
2211  pointMap_(0),
2212  reversePointMap_(0),
2213  pointZone_(0),
2214  retiredPoints_(0),
2215  oldPoints_(0),
2216  faces_(0),
2217  region_(0),
2218  faceOwner_(0),
2219  faceNeighbour_(0),
2220  faceMap_(0),
2221  reverseFaceMap_(0),
2222  faceFromPoint_(0),
2223  faceFromEdge_(0),
2224  flipFaceFlux_(0),
2225  faceZone_(0),
2226  faceZoneFlip_(0),
2227  nActiveFaces_(0),
2228  cellMap_(0),
2229  reverseCellMap_(0),
2230  cellFromPoint_(0),
2231  cellFromEdge_(0),
2232  cellFromFace_(0),
2233  cellZone_(0)
2234 {
2235  addMesh
2236  (
2237  mesh,
2238  identity(mesh.boundaryMesh().size()),
2239  identity(mesh.pointZones().size()),
2240  identity(mesh.faceZones().size()),
2241  identity(mesh.cellZones().size())
2242  );
2243 }
2244 
2245 
2246 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2247 
2249 {
2250  points_.clearStorage();
2251  pointMap_.clearStorage();
2252  reversePointMap_.clearStorage();
2253  pointZone_.clearStorage();
2254  retiredPoints_.clearStorage();
2255  oldPoints_.clearStorage();
2256 
2257  faces_.clearStorage();
2258  region_.clearStorage();
2259  faceOwner_.clearStorage();
2260  faceNeighbour_.clearStorage();
2261  faceMap_.clearStorage();
2262  reverseFaceMap_.clearStorage();
2263  faceFromPoint_.clearStorage();
2264  faceFromEdge_.clearStorage();
2265  flipFaceFlux_.clearStorage();
2266  faceZone_.clearStorage();
2267  faceZoneFlip_.clearStorage();
2268  nActiveFaces_ = 0;
2269 
2270  cellMap_.clearStorage();
2271  reverseCellMap_.clearStorage();
2272  cellZone_.clearStorage();
2273  cellFromPoint_.clearStorage();
2274  cellFromEdge_.clearStorage();
2275  cellFromFace_.clearStorage();
2276 }
2277 
2278 
2281  const polyMesh& mesh,
2282  const labelList& patchMap,
2283  const labelList& pointZoneMap,
2284  const labelList& faceZoneMap,
2285  const labelList& cellZoneMap
2286 )
2287 {
2288  label maxRegion = nPatches_ - 1;
2289  forAll(patchMap, i)
2290  {
2291  maxRegion = max(maxRegion, patchMap[i]);
2292  }
2293  nPatches_ = maxRegion + 1;
2294 
2295 
2296  // Add points
2297  {
2298  const pointField& points = mesh.points();
2299  const pointZoneMesh& pointZones = mesh.pointZones();
2300 
2301  // Extend
2302  points_.setCapacity(points_.size() + points.size());
2303  pointMap_.setCapacity(pointMap_.size() + points.size());
2304  reversePointMap_.setCapacity(reversePointMap_.size() + points.size());
2305  pointZone_.resize(pointZone_.size() + points.size()/100);
2306  // No need to extend oldPoints_
2307 
2308  // Precalc offset zones
2309  labelList newZoneID(points.size(), -1);
2310 
2311  forAll(pointZones, zoneI)
2312  {
2313  const labelList& pointLabels = pointZones[zoneI];
2314 
2315  forAll(pointLabels, j)
2316  {
2317  newZoneID[pointLabels[j]] = pointZoneMap[zoneI];
2318  }
2319  }
2320 
2321  // Add points in mesh order
2322  for (label pointi = 0; pointi < mesh.nPoints(); pointi++)
2323  {
2324  addPoint
2325  (
2326  points[pointi],
2327  pointi,
2328  newZoneID[pointi],
2329  true
2330  );
2331  }
2332  }
2333 
2334  // Add cells
2335  {
2336  const cellZoneMesh& cellZones = mesh.cellZones();
2337 
2338  // Resize
2339 
2340  // Note: polyMesh does not allow retired cells anymore. So allCells
2341  // always equals nCells
2342  label nAllCells = mesh.nCells();
2343 
2344  cellMap_.setCapacity(cellMap_.size() + nAllCells);
2345  reverseCellMap_.setCapacity(reverseCellMap_.size() + nAllCells);
2346  cellFromPoint_.resize(cellFromPoint_.size() + nAllCells/100);
2347  cellFromEdge_.resize(cellFromEdge_.size() + nAllCells/100);
2348  cellFromFace_.resize(cellFromFace_.size() + nAllCells/100);
2349  cellZone_.setCapacity(cellZone_.size() + nAllCells);
2350 
2351 
2352  // Precalc offset zones
2353  labelList newZoneID(nAllCells, -1);
2354 
2355  forAll(cellZones, zoneI)
2356  {
2357  const labelList& cellLabels = cellZones[zoneI];
2358 
2359  forAll(cellLabels, j)
2360  {
2361  label celli = cellLabels[j];
2362 
2363  if (newZoneID[celli] != -1)
2364  {
2366  << "Cell:" << celli
2367  << " centre:" << mesh.cellCentres()[celli]
2368  << " is in two zones:"
2369  << cellZones[newZoneID[celli]].name()
2370  << " and " << cellZones[zoneI].name() << endl
2371  << " This is not supported."
2372  << " Continuing with first zone only." << endl;
2373  }
2374  else
2375  {
2376  newZoneID[celli] = cellZoneMap[zoneI];
2377  }
2378  }
2379  }
2380 
2381  // Add cells in mesh order
2382  for (label celli = 0; celli < nAllCells; celli++)
2383  {
2384  // Add cell from cell
2385  addCell(-1, -1, -1, celli, newZoneID[celli]);
2386  }
2387  }
2388 
2389  // Add faces
2390  {
2391  const polyBoundaryMesh& patches = mesh.boundaryMesh();
2392  const faceList& faces = mesh.faces();
2393  const labelList& faceOwner = mesh.faceOwner();
2394  const labelList& faceNeighbour = mesh.faceNeighbour();
2395  const faceZoneMesh& faceZones = mesh.faceZones();
2396 
2397  // Resize
2398  label nAllFaces = mesh.faces().size();
2399 
2400  faces_.setCapacity(faces_.size() + nAllFaces);
2401  region_.setCapacity(region_.size() + nAllFaces);
2402  faceOwner_.setCapacity(faceOwner_.size() + nAllFaces);
2403  faceNeighbour_.setCapacity(faceNeighbour_.size() + nAllFaces);
2404  faceMap_.setCapacity(faceMap_.size() + nAllFaces);
2405  reverseFaceMap_.setCapacity(reverseFaceMap_.size() + nAllFaces);
2406  faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/100);
2407  faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/100);
2408  flipFaceFlux_.setCapacity(faces_.size() + nAllFaces);
2409  faceZone_.resize(faceZone_.size() + nAllFaces/100);
2410  faceZoneFlip_.setCapacity(faces_.size() + nAllFaces);
2411 
2412 
2413  // Precalc offset zones
2414  labelList newZoneID(nAllFaces, -1);
2415  boolList zoneFlip(nAllFaces, false);
2416 
2417  forAll(faceZones, zoneI)
2418  {
2419  const labelList& faceLabels = faceZones[zoneI];
2420  const boolList& flipMap = faceZones[zoneI].flipMap();
2421 
2422  forAll(faceLabels, j)
2423  {
2424  newZoneID[faceLabels[j]] = faceZoneMap[zoneI];
2425  zoneFlip[faceLabels[j]] = flipMap[j];
2426  }
2427  }
2428 
2429  // Add faces in mesh order
2430 
2431  // 1. Internal faces
2432  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
2433  {
2434  addFace
2435  (
2436  faces[facei],
2437  faceOwner[facei],
2438  faceNeighbour[facei],
2439  -1, // masterPointID
2440  -1, // masterEdgeID
2441  facei, // masterFaceID
2442  false, // flipFaceFlux
2443  -1, // patchID
2444  newZoneID[facei], // zoneID
2445  zoneFlip[facei] // zoneFlip
2446  );
2447  }
2448 
2449  // 2. Patch faces
2450  forAll(patches, patchi)
2451  {
2452  const polyPatch& pp = patches[patchi];
2453 
2454  if (pp.start() != faces_.size())
2455  {
2457  << "Problem : "
2458  << "Patch " << pp.name() << " starts at " << pp.start()
2459  << endl
2460  << "Current face counter at " << faces_.size() << endl
2461  << "Are patches in incremental order?"
2462  << abort(FatalError);
2463  }
2464  forAll(pp, patchFacei)
2465  {
2466  label facei = pp.start() + patchFacei;
2467 
2468  addFace
2469  (
2470  faces[facei],
2471  faceOwner[facei],
2472  -1, // neighbour
2473  -1, // masterPointID
2474  -1, // masterEdgeID
2475  facei, // masterFaceID
2476  false, // flipFaceFlux
2477  patchMap[patchi], // patchID
2478  newZoneID[facei], // zoneID
2479  zoneFlip[facei] // zoneFlip
2480  );
2481  }
2482  }
2483  }
2484 }
2485 
2486 
2489  const label nPoints,
2490  const label nFaces,
2491  const label nCells
2492 )
2493 {
2494  points_.setCapacity(nPoints);
2495  pointMap_.setCapacity(nPoints);
2496  reversePointMap_.setCapacity(nPoints);
2497  pointZone_.resize(pointZone_.size() + nPoints/100);
2498 
2499  faces_.setCapacity(nFaces);
2500  region_.setCapacity(nFaces);
2501  faceOwner_.setCapacity(nFaces);
2502  faceNeighbour_.setCapacity(nFaces);
2503  faceMap_.setCapacity(nFaces);
2504  reverseFaceMap_.setCapacity(nFaces);
2505  faceFromPoint_.resize(faceFromPoint_.size() + nFaces/100);
2506  faceFromEdge_.resize(faceFromEdge_.size() + nFaces/100);
2507  flipFaceFlux_.setCapacity(nFaces);
2508  faceZone_.resize(faceZone_.size() + nFaces/100);
2509  faceZoneFlip_.setCapacity(nFaces);
2510 
2511  cellMap_.setCapacity(nCells);
2512  reverseCellMap_.setCapacity(nCells);
2513  cellFromPoint_.resize(cellFromPoint_.size() + nCells/100);
2514  cellFromEdge_.resize(cellFromEdge_.size() + nCells/100);
2515  cellFromFace_.resize(cellFromFace_.size() + nCells/100);
2516  cellZone_.setCapacity(nCells);
2517 }
2518 
2519 
2521 {
2522  if (isType<polyAddPoint>(action))
2523  {
2524  const polyAddPoint& pap = refCast<const polyAddPoint>(action);
2525 
2526  return addPoint
2527  (
2528  pap.newPoint(),
2529  pap.masterPointID(),
2530  pap.zoneID(),
2531  pap.inCell()
2532  );
2533  }
2534  else if (isType<polyModifyPoint>(action))
2535  {
2536  const polyModifyPoint& pmp = refCast<const polyModifyPoint>(action);
2537 
2538  modifyPoint
2539  (
2540  pmp.pointID(),
2541  pmp.newPoint(),
2542  pmp.zoneID(),
2543  pmp.inCell()
2544  );
2545 
2546  return -1;
2547  }
2548  else if (isType<polyRemovePoint>(action))
2549  {
2550  const polyRemovePoint& prp = refCast<const polyRemovePoint>(action);
2551 
2552  removePoint(prp.pointID(), prp.mergePointID());
2553 
2554  return -1;
2555  }
2556  else if (isType<polyAddFace>(action))
2557  {
2558  const polyAddFace& paf = refCast<const polyAddFace>(action);
2559 
2560  return addFace
2561  (
2562  paf.newFace(),
2563  paf.owner(),
2564  paf.neighbour(),
2565  paf.masterPointID(),
2566  paf.masterEdgeID(),
2567  paf.masterFaceID(),
2568  paf.flipFaceFlux(),
2569  paf.patchID(),
2570  paf.zoneID(),
2571  paf.zoneFlip()
2572  );
2573  }
2574  else if (isType<polyModifyFace>(action))
2575  {
2576  const polyModifyFace& pmf = refCast<const polyModifyFace>(action);
2577 
2578  modifyFace
2579  (
2580  pmf.newFace(),
2581  pmf.faceID(),
2582  pmf.owner(),
2583  pmf.neighbour(),
2584  pmf.flipFaceFlux(),
2585  pmf.patchID(),
2586  pmf.zoneID(),
2587  pmf.zoneFlip()
2588  );
2589 
2590  return -1;
2591  }
2592  else if (isType<polyRemoveFace>(action))
2593  {
2594  const polyRemoveFace& prf = refCast<const polyRemoveFace>(action);
2595 
2596  removeFace(prf.faceID(), prf.mergeFaceID());
2597 
2598  return -1;
2599  }
2600  else if (isType<polyAddCell>(action))
2601  {
2602  const polyAddCell& pac = refCast<const polyAddCell>(action);
2603 
2604  return addCell
2605  (
2606  pac.masterPointID(),
2607  pac.masterEdgeID(),
2608  pac.masterFaceID(),
2609  pac.masterCellID(),
2610  pac.zoneID()
2611  );
2612  }
2613  else if (isType<polyModifyCell>(action))
2614  {
2615  const polyModifyCell& pmc = refCast<const polyModifyCell>(action);
2616 
2617  if (pmc.removeFromZone())
2618  {
2619  modifyCell(pmc.cellID(), -1);
2620  }
2621  else
2622  {
2623  modifyCell(pmc.cellID(), pmc.zoneID());
2624  }
2625 
2626  return -1;
2627  }
2628  else if (isType<polyRemoveCell>(action))
2629  {
2630  const polyRemoveCell& prc = refCast<const polyRemoveCell>(action);
2631 
2632  removeCell(prc.cellID(), prc.mergeCellID());
2633 
2634  return -1;
2635  }
2636  else
2637  {
2639  << "Unknown type of topoChange: " << action.type()
2640  << abort(FatalError);
2641 
2642  // Dummy return to keep compiler happy
2643  return -1;
2644  }
2645 }
2646 
2647 
2650  const point& pt,
2651  const label masterPointID,
2652  const label zoneID,
2653  const bool inCell
2654 )
2655 {
2656  label pointi = points_.size();
2657 
2658  points_.append(pt);
2659  pointMap_.append(masterPointID);
2660  reversePointMap_.append(pointi);
2661 
2662  if (zoneID >= 0)
2663  {
2664  pointZone_.insert(pointi, zoneID);
2665  }
2666 
2667  if (!inCell)
2668  {
2669  retiredPoints_.insert(pointi);
2670  }
2671 
2672  return pointi;
2673 }
2674 
2675 
2678  const label pointi,
2679  const point& pt,
2680  const label newZoneID,
2681  const bool inCell
2682 )
2683 {
2684  if (pointi < 0 || pointi >= points_.size())
2685  {
2687  << "illegal point label " << pointi << endl
2688  << "Valid point labels are 0 .. " << points_.size()-1
2689  << abort(FatalError);
2690  }
2691  if (pointRemoved(pointi) || pointMap_[pointi] == -1)
2692  {
2694  << "point " << pointi << " already marked for removal"
2695  << abort(FatalError);
2696  }
2697  points_[pointi] = pt;
2698 
2699  Map<label>::iterator pointFnd = pointZone_.find(pointi);
2700 
2701  if (pointFnd != pointZone_.end())
2702  {
2703  if (newZoneID >= 0)
2704  {
2705  pointFnd() = newZoneID;
2706  }
2707  else
2708  {
2709  pointZone_.erase(pointFnd);
2710  }
2711  }
2712  else if (newZoneID >= 0)
2713  {
2714  pointZone_.insert(pointi, newZoneID);
2715  }
2716 
2717  if (inCell)
2718  {
2719  retiredPoints_.erase(pointi);
2720  }
2721  else
2722  {
2723  retiredPoints_.insert(pointi);
2724  }
2725 
2726  oldPoints_.erase(pointi);
2727 }
2728 
2729 
2732  const point& pt,
2733  const point& oldPt,
2734  const label masterPointID,
2735  const label zoneID
2736 )
2737 {
2738  label pointi = points_.size();
2739 
2740  points_.append(pt);
2741  pointMap_.append(masterPointID);
2742  reversePointMap_.append(pointi);
2743 
2744  if (zoneID >= 0)
2745  {
2746  pointZone_.insert(pointi, zoneID);
2747  }
2748 
2749  oldPoints_.insert(pointi, oldPt);
2750 
2751  return pointi;
2752 }
2753 
2754 
2757  const label pointi,
2758  const point& pt,
2759  const point& oldPt,
2760  const label newZoneID
2761 )
2762 {
2763  if (pointi < 0 || pointi >= points_.size())
2764  {
2766  << "illegal point label " << pointi << endl
2767  << "Valid point labels are 0 .. " << points_.size()-1
2768  << abort(FatalError);
2769  }
2770  if (pointRemoved(pointi) || pointMap_[pointi] == -1)
2771  {
2773  << "point " << pointi << " already marked for removal"
2774  << abort(FatalError);
2775  }
2776  points_[pointi] = pt;
2777 
2778  Map<label>::iterator pointFnd = pointZone_.find(pointi);
2779 
2780  if (pointFnd != pointZone_.end())
2781  {
2782  if (newZoneID >= 0)
2783  {
2784  pointFnd() = newZoneID;
2785  }
2786  else
2787  {
2788  pointZone_.erase(pointFnd);
2789  }
2790  }
2791  else if (newZoneID >= 0)
2792  {
2793  pointZone_.insert(pointi, newZoneID);
2794  }
2795 
2796  // Always active
2797  retiredPoints_.erase(pointi);
2798 
2799  // Always provided old point
2800  oldPoints_.set(pointi, oldPt);
2801 }
2802 
2803 
2805 {
2806  if (newPoints.size() != points_.size())
2807  {
2809  << "illegal pointField size." << endl
2810  << "Size:" << newPoints.size() << endl
2811  << "Points in mesh:" << points_.size()
2812  << abort(FatalError);
2813  }
2814 
2815  forAll(points_, pointi)
2816  {
2817  points_[pointi] = newPoints[pointi];
2818  }
2819 }
2820 
2821 
2824  const label pointi,
2825  const label mergePointi
2826 )
2827 {
2828  if (pointi < 0 || pointi >= points_.size())
2829  {
2831  << "illegal point label " << pointi << endl
2832  << "Valid point labels are 0 .. " << points_.size()-1
2833  << abort(FatalError);
2834  }
2835 
2836  if
2837  (
2838  strict_
2839  && (pointRemoved(pointi) || pointMap_[pointi] == -1)
2840  )
2841  {
2843  << "point " << pointi << " already marked for removal" << nl
2844  << "Point:" << points_[pointi] << " pointMap:" << pointMap_[pointi]
2845  << abort(FatalError);
2846  }
2847 
2848  if (pointi == mergePointi)
2849  {
2851  << "Cannot remove/merge point " << pointi << " onto itself."
2852  << abort(FatalError);
2853  }
2854 
2855  points_[pointi] = point::max;
2856  pointMap_[pointi] = -1;
2857  if (mergePointi >= 0)
2858  {
2859  reversePointMap_[pointi] = -mergePointi-2;
2860  }
2861  else
2862  {
2863  reversePointMap_[pointi] = -1;
2864  }
2865  pointZone_.erase(pointi);
2866  retiredPoints_.erase(pointi);
2867  oldPoints_.erase(pointi);
2868 }
2869 
2870 
2873  const face& f,
2874  const label own,
2875  const label nei,
2876  const label masterPointID,
2877  const label masterEdgeID,
2878  const label masterFaceID,
2879  const bool flipFaceFlux,
2880  const label patchID,
2881  const label zoneID,
2882  const bool zoneFlip
2883 )
2884 {
2885  // Check validity
2886  if (debug)
2887  {
2888  checkFace(f, -1, own, nei, patchID, zoneID);
2889  }
2890 
2891  label facei = faces_.size();
2892 
2893  faces_.append(f);
2894  region_.append(patchID);
2895  faceOwner_.append(own);
2896  faceNeighbour_.append(nei);
2897 
2898  if (masterPointID >= 0)
2899  {
2900  faceMap_.append(-1);
2901  faceFromPoint_.insert(facei, masterPointID);
2902  }
2903  else if (masterEdgeID >= 0)
2904  {
2905  faceMap_.append(-1);
2906  faceFromEdge_.insert(facei, masterEdgeID);
2907  }
2908  else if (masterFaceID >= 0)
2909  {
2910  faceMap_.append(masterFaceID);
2911  }
2912  else
2913  {
2914  // Allow inflate-from-nothing?
2915  // FatalErrorInFunction
2916  // << "Need to specify a master point, edge or face"
2917  // << "face:" << f << " own:" << own << " nei:" << nei
2918  // << abort(FatalError);
2919  faceMap_.append(-1);
2920  }
2921  reverseFaceMap_.append(facei);
2922 
2923  flipFaceFlux_[facei] = (flipFaceFlux ? 1 : 0);
2924 
2925  if (zoneID >= 0)
2926  {
2927  faceZone_.insert(facei, zoneID);
2928  }
2929  faceZoneFlip_[facei] = (zoneFlip ? 1 : 0);
2930 
2931  return facei;
2932 }
2933 
2934 
2937  const face& f,
2938  const label facei,
2939  const label own,
2940  const label nei,
2941  const bool flipFaceFlux,
2942  const label patchID,
2943  const label zoneID,
2944  const bool zoneFlip
2945 )
2946 {
2947  // Check validity
2948  if (debug)
2949  {
2950  checkFace(f, facei, own, nei, patchID, zoneID);
2951  }
2952 
2953  faces_[facei] = f;
2954  faceOwner_[facei] = own;
2955  faceNeighbour_[facei] = nei;
2956  region_[facei] = patchID;
2957 
2958  flipFaceFlux_[facei] = (flipFaceFlux ? 1 : 0);
2959 
2960  Map<label>::iterator faceFnd = faceZone_.find(facei);
2961 
2962  if (faceFnd != faceZone_.end())
2963  {
2964  if (zoneID >= 0)
2965  {
2966  faceFnd() = zoneID;
2967  }
2968  else
2969  {
2970  faceZone_.erase(faceFnd);
2971  }
2972  }
2973  else if (zoneID >= 0)
2974  {
2975  faceZone_.insert(facei, zoneID);
2976  }
2977  faceZoneFlip_[facei] = (zoneFlip ? 1 : 0);
2978 }
2979 
2980 
2981 void Foam::polyTopoChange::removeFace(const label facei, const label mergeFacei)
2982 {
2983  if (facei < 0 || facei >= faces_.size())
2984  {
2986  << "illegal face label " << facei << endl
2987  << "Valid face labels are 0 .. " << faces_.size()-1
2988  << abort(FatalError);
2989  }
2990 
2991  if
2992  (
2993  strict_
2994  && (faceRemoved(facei) || faceMap_[facei] == -1)
2995  )
2996  {
2998  << "face " << facei
2999  << " already marked for removal"
3000  << abort(FatalError);
3001  }
3002 
3003  faces_[facei].setSize(0);
3004  region_[facei] = -1;
3005  faceOwner_[facei] = -1;
3006  faceNeighbour_[facei] = -1;
3007  faceMap_[facei] = -1;
3008  if (mergeFacei >= 0)
3009  {
3010  reverseFaceMap_[facei] = -mergeFacei-2;
3011  }
3012  else
3013  {
3014  reverseFaceMap_[facei] = -1;
3015  }
3016  faceFromEdge_.erase(facei);
3017  faceFromPoint_.erase(facei);
3018  flipFaceFlux_[facei] = 0;
3019  faceZone_.erase(facei);
3020  faceZoneFlip_[facei] = 0;
3021 }
3022 
3023 
3026  const label masterPointID,
3027  const label masterEdgeID,
3028  const label masterFaceID,
3029  const label masterCellID,
3030  const label zoneID
3031 )
3032 {
3033  label celli = cellMap_.size();
3034 
3035  if (masterPointID >= 0)
3036  {
3037  cellMap_.append(-1);
3038  cellFromPoint_.insert(celli, masterPointID);
3039  }
3040  else if (masterEdgeID >= 0)
3041  {
3042  cellMap_.append(-1);
3043  cellFromEdge_.insert(celli, masterEdgeID);
3044  }
3045  else if (masterFaceID >= 0)
3046  {
3047  cellMap_.append(-1);
3048  cellFromFace_.insert(celli, masterFaceID);
3049  }
3050  else
3051  {
3052  cellMap_.append(masterCellID);
3053  }
3054  reverseCellMap_.append(celli);
3055  cellZone_.append(zoneID);
3056 
3057  return celli;
3058 }
3059 
3060 
3063  const label celli,
3064  const label zoneID
3065 )
3066 {
3067  cellZone_[celli] = zoneID;
3068 }
3069 
3070 
3071 void Foam::polyTopoChange::removeCell(const label celli, const label mergeCelli)
3072 {
3073  if (celli < 0 || celli >= cellMap_.size())
3074  {
3076  << "illegal cell label " << celli << endl
3077  << "Valid cell labels are 0 .. " << cellMap_.size()-1
3078  << abort(FatalError);
3079  }
3080 
3081  if (strict_ && cellMap_[celli] == -2)
3082  {
3084  << "cell " << celli
3085  << " already marked for removal"
3086  << abort(FatalError);
3087  }
3088 
3089  cellMap_[celli] = -2;
3090  if (mergeCelli >= 0)
3091  {
3092  reverseCellMap_[celli] = -mergeCelli-2;
3093  }
3094  else
3095  {
3096  reverseCellMap_[celli] = -1;
3097  }
3098  cellFromPoint_.erase(celli);
3099  cellFromEdge_.erase(celli);
3100  cellFromFace_.erase(celli);
3101  cellZone_[celli] = -1;
3102 }
3103 
3104 
3107  polyMesh& mesh,
3108  const bool inflate,
3109  const bool syncParallel,
3110  const bool orderCells,
3111  const bool orderPoints
3112 )
3113 {
3114  if (debug)
3115  {
3116  Pout<< "polyTopoChange::changeMesh"
3117  << "(polyMesh&, const bool, const bool, const bool, const bool)"
3118  << endl;
3119  }
3120 
3121  if (debug)
3122  {
3123  Pout<< "Old mesh:" << nl;
3124  writeMeshStats(mesh, Pout);
3125  }
3126 
3127  // new mesh points
3128  pointField newPoints;
3129  // number of internal points
3130  label nInternalPoints;
3131  // patch slicing
3132  labelList patchSizes;
3133  labelList patchStarts;
3134  // inflate maps
3135  List<objectMap> pointsFromPoints;
3136  List<objectMap> facesFromPoints;
3137  List<objectMap> facesFromEdges;
3138  List<objectMap> facesFromFaces;
3139  List<objectMap> cellsFromPoints;
3140  List<objectMap> cellsFromEdges;
3141  List<objectMap> cellsFromFaces;
3142  List<objectMap> cellsFromCells;
3143  // old mesh info
3144  List<Map<label>> oldPatchMeshPointMaps;
3145  labelList oldPatchNMeshPoints;
3146  labelList oldPatchStarts;
3147  List<Map<label>> oldFaceZoneMeshPointMaps;
3148 
3149  // Compact, reorder patch faces and calculate mesh/patch maps.
3150  compactAndReorder
3151  (
3152  mesh,
3153  syncParallel,
3154  orderCells,
3155  orderPoints,
3156 
3157  nInternalPoints,
3158  newPoints,
3159  patchSizes,
3160  patchStarts,
3161  pointsFromPoints,
3162  facesFromPoints,
3163  facesFromEdges,
3164  facesFromFaces,
3165  cellsFromPoints,
3166  cellsFromEdges,
3167  cellsFromFaces,
3168  cellsFromCells,
3169  oldPatchMeshPointMaps,
3170  oldPatchNMeshPoints,
3171  oldPatchStarts,
3172  oldFaceZoneMeshPointMaps
3173  );
3174 
3175  const label nOldPoints(mesh.nPoints());
3176  const label nOldFaces(mesh.nFaces());
3177  const label nOldCells(mesh.nCells());
3178  autoPtr<scalarField> oldCellVolumes(new scalarField(mesh.cellVolumes()));
3179 
3180 
3181  // Change the mesh
3182  // ~~~~~~~~~~~~~~~
3183  // This will invalidate any addressing so better make sure you have
3184  // all the information you need!!!
3185 
3186  if (inflate)
3187  {
3188  // Keep (renumbered) mesh points, store new points in map for inflation
3189  // (appended points (i.e. from nowhere) get value as provided in
3190  // addPoints)
3191  pointField renumberedMeshPoints(newPoints.size());
3192 
3193  forAll(pointMap_, newPointi)
3194  {
3195  Map<point>::const_iterator iter = oldPoints_.find(newPointi);
3196  if (iter != oldPoints_.end())
3197  {
3198  renumberedMeshPoints[newPointi] = iter();
3199  }
3200  else
3201  {
3202  label oldPointi = pointMap_[newPointi];
3203 
3204  if (oldPointi >= 0)
3205  {
3206  renumberedMeshPoints[newPointi] = mesh.points()[oldPointi];
3207  }
3208  else
3209  {
3210  renumberedMeshPoints[newPointi] = vector::zero;
3211  }
3212  }
3213  }
3214 
3215  mesh.resetPrimitives
3216  (
3217  move(renumberedMeshPoints),
3218  move(faces_),
3219  move(faceOwner_),
3220  move(faceNeighbour_),
3221  patchSizes,
3222  patchStarts,
3223  syncParallel
3224  );
3225 
3226  mesh.topoChanging(true);
3227  // Note: could already set moving flag as well
3228  // mesh.moving(true);
3229  }
3230  else
3231  {
3232  // Set new points.
3233  mesh.resetPrimitives
3234  (
3235  move(newPoints),
3236  move(faces_),
3237  move(faceOwner_),
3238  move(faceNeighbour_),
3239  patchSizes,
3240  patchStarts,
3241  syncParallel
3242  );
3243  mesh.topoChanging(true);
3244  }
3245 
3246  // Clear out primitives
3247  {
3248  retiredPoints_.clearStorage();
3249  oldPoints_.clearStorage();
3250  region_.clearStorage();
3251  }
3252 
3253 
3254  if (debug)
3255  {
3256  // Some stats on changes
3257  label nAdd, nInflate, nMerge, nRemove;
3258  countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3259  Pout<< "Points:"
3260  << " added(from point):" << nAdd
3261  << " added(from nothing):" << nInflate
3262  << " merged(into other point):" << nMerge
3263  << " removed:" << nRemove
3264  << nl;
3265 
3266  countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3267  Pout<< "Faces:"
3268  << " added(from face):" << nAdd
3269  << " added(inflated):" << nInflate
3270  << " merged(into other face):" << nMerge
3271  << " removed:" << nRemove
3272  << nl;
3273 
3274  countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3275  Pout<< "Cells:"
3276  << " added(from cell):" << nAdd
3277  << " added(inflated):" << nInflate
3278  << " merged(into other cell):" << nMerge
3279  << " removed:" << nRemove
3280  << nl
3281  << endl;
3282  }
3283 
3284  if (debug)
3285  {
3286  Pout<< "New mesh:" << nl;
3287  writeMeshStats(mesh, Pout);
3288  }
3289 
3290 
3291  // Zones
3292  // ~~~~~
3293 
3294  // Inverse of point/face/cell zone addressing.
3295  // For every preserved point/face/cells in zone give the old position.
3296  // For added points, the index is set to -1
3297  labelListList pointZoneMap(mesh.pointZones().size());
3298  labelListList faceZoneFaceMap(mesh.faceZones().size());
3299  labelListList cellZoneMap(mesh.cellZones().size());
3300 
3301  resetZones(mesh, mesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3302 
3303  // Clear zone info
3304  {
3305  pointZone_.clearStorage();
3306  faceZone_.clearStorage();
3307  faceZoneFlip_.clearStorage();
3308  cellZone_.clearStorage();
3309  }
3310 
3311 
3312  // Patch point renumbering
3313  // For every preserved point on a patch give the old position.
3314  // For added points, the index is set to -1
3315  labelListList patchPointMap(mesh.boundaryMesh().size());
3316  calcPatchPointMap
3317  (
3318  oldPatchMeshPointMaps,
3319  mesh.boundaryMesh(),
3320  patchPointMap
3321  );
3322 
3323  // Create the face zone mesh point renumbering
3324  labelListList faceZonePointMap(mesh.faceZones().size());
3325  calcFaceZonePointMap(mesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3326 
3327  labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3328 
3329  return autoPtr<mapPolyMesh>
3330  (
3331  new mapPolyMesh
3332  (
3333  mesh,
3334  nOldPoints,
3335  nOldFaces,
3336  nOldCells,
3337 
3338  pointMap_,
3339  pointsFromPoints,
3340 
3341  faceMap_,
3342  facesFromPoints,
3343  facesFromEdges,
3344  facesFromFaces,
3345 
3346  cellMap_,
3347  cellsFromPoints,
3348  cellsFromEdges,
3349  cellsFromFaces,
3350  cellsFromCells,
3351 
3352  reversePointMap_,
3353  reverseFaceMap_,
3354  reverseCellMap_,
3355 
3356  flipFaceFluxSet,
3357 
3358  patchPointMap,
3359 
3360  pointZoneMap,
3361 
3362  faceZonePointMap,
3363  faceZoneFaceMap,
3364  cellZoneMap,
3365 
3366  newPoints, // if empty signals no inflation.
3367  oldPatchStarts,
3368  oldPatchNMeshPoints,
3369 
3370  oldCellVolumes,
3371 
3372  true // steal storage.
3373  )
3374  );
3375 
3376  // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3377  // be invalid.
3378 }
3379 
3380 
3383  autoPtr<fvMesh>& newMeshPtr,
3384  const IOobject& io,
3385  const polyMesh& mesh,
3386  const bool syncParallel,
3387  const bool orderCells,
3388  const bool orderPoints
3389 )
3390 {
3391  if (debug)
3392  {
3393  Pout<< "polyTopoChange::changeMesh"
3394  << "(autoPtr<fvMesh>&, const IOobject&, const fvMesh&"
3395  << ", const bool, const bool, const bool)"
3396  << endl;
3397  }
3398 
3399  if (debug)
3400  {
3401  Pout<< "Old mesh:" << nl;
3402  writeMeshStats(mesh, Pout);
3403  }
3404 
3405  // new mesh points
3406  pointField newPoints;
3407  // number of internal points
3408  label nInternalPoints;
3409  // patch slicing
3410  labelList patchSizes;
3411  labelList patchStarts;
3412  // inflate maps
3413  List<objectMap> pointsFromPoints;
3414  List<objectMap> facesFromPoints;
3415  List<objectMap> facesFromEdges;
3416  List<objectMap> facesFromFaces;
3417  List<objectMap> cellsFromPoints;
3418  List<objectMap> cellsFromEdges;
3419  List<objectMap> cellsFromFaces;
3420  List<objectMap> cellsFromCells;
3421 
3422  // old mesh info
3423  List<Map<label>> oldPatchMeshPointMaps;
3424  labelList oldPatchNMeshPoints;
3425  labelList oldPatchStarts;
3426  List<Map<label>> oldFaceZoneMeshPointMaps;
3427 
3428  // Compact, reorder patch faces and calculate mesh/patch maps.
3429  compactAndReorder
3430  (
3431  mesh,
3432  syncParallel,
3433  orderCells,
3434  orderPoints,
3435 
3436  nInternalPoints,
3437  newPoints,
3438  patchSizes,
3439  patchStarts,
3440  pointsFromPoints,
3441  facesFromPoints,
3442  facesFromEdges,
3443  facesFromFaces,
3444  cellsFromPoints,
3445  cellsFromEdges,
3446  cellsFromFaces,
3447  cellsFromCells,
3448  oldPatchMeshPointMaps,
3449  oldPatchNMeshPoints,
3450  oldPatchStarts,
3451  oldFaceZoneMeshPointMaps
3452  );
3453 
3454  const label nOldPoints(mesh.nPoints());
3455  const label nOldFaces(mesh.nFaces());
3456  const label nOldCells(mesh.nCells());
3457  autoPtr<scalarField> oldCellVolumes(new scalarField(mesh.cellVolumes()));
3458 
3459 
3460  // Create the mesh
3461  // ~~~~~~~~~~~~~~~
3462 
3463  IOobject noReadIO(io);
3464  noReadIO.readOpt() = IOobject::NO_READ;
3465  newMeshPtr.reset
3466  (
3467  new fvMesh
3468  (
3469  noReadIO,
3470  move(newPoints),
3471  move(faces_),
3472  move(faceOwner_),
3473  move(faceNeighbour_)
3474  )
3475  );
3476  fvMesh& newMesh = newMeshPtr();
3477 
3478  // Clear out primitives
3479  {
3480  retiredPoints_.clearStorage();
3481  oldPoints_.clearStorage();
3482  region_.clearStorage();
3483  }
3484 
3485 
3486  if (debug)
3487  {
3488  // Some stats on changes
3489  label nAdd, nInflate, nMerge, nRemove;
3490  countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3491  Pout<< "Points:"
3492  << " added(from point):" << nAdd
3493  << " added(from nothing):" << nInflate
3494  << " merged(into other point):" << nMerge
3495  << " removed:" << nRemove
3496  << nl;
3497 
3498  countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3499  Pout<< "Faces:"
3500  << " added(from face):" << nAdd
3501  << " added(inflated):" << nInflate
3502  << " merged(into other face):" << nMerge
3503  << " removed:" << nRemove
3504  << nl;
3505 
3506  countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3507  Pout<< "Cells:"
3508  << " added(from cell):" << nAdd
3509  << " added(inflated):" << nInflate
3510  << " merged(into other cell):" << nMerge
3511  << " removed:" << nRemove
3512  << nl
3513  << endl;
3514  }
3515 
3516 
3517  {
3518  const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
3519 
3520  List<polyPatch*> newBoundary(oldPatches.size());
3521 
3522  forAll(oldPatches, patchi)
3523  {
3524  newBoundary[patchi] = oldPatches[patchi].clone
3525  (
3526  newMesh.boundaryMesh(),
3527  patchi,
3528  patchSizes[patchi],
3529  patchStarts[patchi]
3530  ).ptr();
3531  }
3532  newMesh.addFvPatches(newBoundary);
3533  }
3534 
3535 
3536  // Zones
3537  // ~~~~~
3538 
3539  // Start off from empty zones.
3540  const pointZoneMesh& oldPointZones = mesh.pointZones();
3541  List<pointZone*> pZonePtrs(oldPointZones.size());
3542  {
3543  forAll(oldPointZones, i)
3544  {
3545  pZonePtrs[i] = new pointZone
3546  (
3547  oldPointZones[i].name(),
3548  labelList(0),
3549  i,
3550  newMesh.pointZones()
3551  );
3552  }
3553  }
3554 
3555  const faceZoneMesh& oldFaceZones = mesh.faceZones();
3556  List<faceZone*> fZonePtrs(oldFaceZones.size());
3557  {
3558  forAll(oldFaceZones, i)
3559  {
3560  fZonePtrs[i] = new faceZone
3561  (
3562  oldFaceZones[i].name(),
3563  labelList(0),
3564  boolList(0),
3565  i,
3566  newMesh.faceZones()
3567  );
3568  }
3569  }
3570 
3571  const cellZoneMesh& oldCellZones = mesh.cellZones();
3572  List<cellZone*> cZonePtrs(oldCellZones.size());
3573  {
3574  forAll(oldCellZones, i)
3575  {
3576  cZonePtrs[i] = new cellZone
3577  (
3578  oldCellZones[i].name(),
3579  labelList(0),
3580  i,
3581  newMesh.cellZones()
3582  );
3583  }
3584  }
3585 
3586  newMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
3587 
3588  // Inverse of point/face/cell zone addressing.
3589  // For every preserved point/face/cells in zone give the old position.
3590  // For added points, the index is set to -1
3591  labelListList pointZoneMap(mesh.pointZones().size());
3592  labelListList faceZoneFaceMap(mesh.faceZones().size());
3593  labelListList cellZoneMap(mesh.cellZones().size());
3594 
3595  resetZones(mesh, newMesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3596 
3597  // Clear zone info
3598  {
3599  pointZone_.clearStorage();
3600  faceZone_.clearStorage();
3601  faceZoneFlip_.clearStorage();
3602  cellZone_.clearStorage();
3603  }
3604 
3605  // Patch point renumbering
3606  // For every preserved point on a patch give the old position.
3607  // For added points, the index is set to -1
3608  labelListList patchPointMap(newMesh.boundaryMesh().size());
3609  calcPatchPointMap
3610  (
3611  oldPatchMeshPointMaps,
3612  newMesh.boundaryMesh(),
3613  patchPointMap
3614  );
3615 
3616  // Create the face zone mesh point renumbering
3617  labelListList faceZonePointMap(newMesh.faceZones().size());
3618  calcFaceZonePointMap(newMesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3619 
3620  if (debug)
3621  {
3622  Pout<< "New mesh:" << nl;
3623  writeMeshStats(mesh, Pout);
3624  }
3625 
3626  labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3627 
3628  return autoPtr<mapPolyMesh>
3629  (
3630  new mapPolyMesh
3631  (
3632  newMesh,
3633  nOldPoints,
3634  nOldFaces,
3635  nOldCells,
3636 
3637  pointMap_,
3638  pointsFromPoints,
3639 
3640  faceMap_,
3641  facesFromPoints,
3642  facesFromEdges,
3643  facesFromFaces,
3644 
3645  cellMap_,
3646  cellsFromPoints,
3647  cellsFromEdges,
3648  cellsFromFaces,
3649  cellsFromCells,
3650 
3651  reversePointMap_,
3652  reverseFaceMap_,
3653  reverseCellMap_,
3654 
3655  flipFaceFluxSet,
3656 
3657  patchPointMap,
3658 
3659  pointZoneMap,
3660 
3661  faceZonePointMap,
3662  faceZoneFaceMap,
3663  cellZoneMap,
3664 
3665  newPoints, // if empty signals no inflation.
3666  oldPatchStarts,
3667  oldPatchNMeshPoints,
3668  oldCellVolumes,
3669  true // steal storage.
3670  )
3671  );
3672 
3673  // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3674  // be invalid.
3675 }
3676 
3677 
3678 // ************************************************************************* //
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:268
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
#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:295
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
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:476
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)
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:319
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:163
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:256
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:108
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
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:208
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
const pointZoneMesh & pointZones() const
Return point zone mesh.
Definition: polyMesh.H:470
label nPoints
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const cellZoneMesh & cellZones() const
Return cell zone mesh.
Definition: polyMesh.H:482
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.
ZoneMesh< pointZone, polyMesh > pointZoneMesh
A ZoneMesh with the type pointZone.
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1156
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
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:265
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.
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 neighour 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
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:303
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:345
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.
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.
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