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