FreeFOAM The Cross-Platform CFD Toolkit
combineFaces.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "combineFaces.H"
27 #include <OpenFOAM/polyMesh.H>
28 #include "polyTopoChange.H"
34 #include <OpenFOAM/syncTools.H>
35 #include <meshTools/meshTools.H>
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41 
42 defineTypeNameAndDebug(combineFaces, 0);
43 
44 }
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 // Test face for (almost) convexeness. Allows certain convexity before
49 // complaining.
50 // For every two consecutive edges calculate the normal. If it differs too
51 // much from the average face normal complain.
52 bool Foam::combineFaces::convexFace
53 (
54  const scalar minConcaveCos,
55  const pointField& points,
56  const face& f
57 )
58 {
59  // Get outwards pointing normal of f.
60  vector n = f.normal(points);
61  n /= mag(n);
62 
63  // Get edge from f[0] to f[size-1];
64  vector ePrev(points[f[0]] - points[f[f.size()-1]]);
65  scalar magEPrev = mag(ePrev);
66  ePrev /= magEPrev + VSMALL;
67 
68  forAll(f, fp0)
69  {
70  // Get vertex after fp
71  label fp1 = f.fcIndex(fp0);
72 
73  // Normalized vector between two consecutive points
74  vector e10(points[f[fp1]] - points[f[fp0]]);
75  scalar magE10 = mag(e10);
76  e10 /= magE10 + VSMALL;
77 
78  if (magEPrev > SMALL && magE10 > SMALL)
79  {
80  vector edgeNormal = ePrev ^ e10;
81 
82  if ((edgeNormal & n) < 0)
83  {
84  // Concave. Check angle.
85  if ((ePrev & e10) < minConcaveCos)
86  {
87  return false;
88  }
89  }
90  }
91 
92  ePrev = e10;
93  magEPrev = magE10;
94  }
95 
96  // Not a single internal angle is concave so face is convex.
97  return true;
98 }
99 
100 
101 // Determines if set of faces is valid to collapse into single face.
102 bool Foam::combineFaces::validFace
103 (
104  const scalar minConcaveCos,
105  const indirectPrimitivePatch& bigFace
106 )
107 {
108  // Get outside vertices (in local vertex numbering)
109  const labelListList& edgeLoops = bigFace.edgeLoops();
110 
111  if (edgeLoops.size() > 1)
112  {
113  return false;
114  }
115 
116  // Check for convexness
117  face f(getOutsideFace(bigFace));
118 
119  return convexFace(minConcaveCos, bigFace.points(), f);
120 }
121 
122 
123 void Foam::combineFaces::regioniseFaces
124 (
125  const scalar minCos,
126  const label cellI,
127  const labelList& cEdges,
128  Map<label>& faceRegion
129 ) const
130 {
131  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
132 
133  forAll(cEdges, i)
134  {
135  label edgeI = cEdges[i];
136 
137  label f0, f1;
138  meshTools::getEdgeFaces(mesh_, cellI, edgeI, f0, f1);
139 
140  label p0 = patches.whichPatch(f0);
141  label p1 = patches.whichPatch(f1);
142 
143  // Face can be merged if
144  // - same non-coupled patch
145  // - small angle
146  if (p0 != -1 && p0 == p1 && !patches[p0].coupled())
147  {
148  vector f0Normal = mesh_.faceAreas()[f0];
149  f0Normal /= mag(f0Normal);
150  vector f1Normal = mesh_.faceAreas()[f1];
151  f1Normal /= mag(f1Normal);
152 
153  if ((f0Normal&f1Normal) > minCos)
154  {
155  Map<label>::const_iterator f0Fnd = faceRegion.find(f0);
156 
157  label region0 = -1;
158  if (f0Fnd != faceRegion.end())
159  {
160  region0 = f0Fnd();
161  }
162 
163  Map<label>::const_iterator f1Fnd = faceRegion.find(f1);
164 
165  label region1 = -1;
166  if (f1Fnd != faceRegion.end())
167  {
168  region1 = f1Fnd();
169  }
170 
171  if (region0 == -1)
172  {
173  if (region1 == -1)
174  {
175  label useRegion = faceRegion.size();
176  faceRegion.insert(f0, useRegion);
177  faceRegion.insert(f1, useRegion);
178  }
179  else
180  {
181  faceRegion.insert(f0, region1);
182  }
183  }
184  else
185  {
186  if (region1 == -1)
187  {
188  faceRegion.insert(f1, region0);
189  }
190  else if (region0 != region1)
191  {
192  // Merge the two regions
193  label useRegion = min(region0, region1);
194  label freeRegion = max(region0, region1);
195 
196  forAllIter(Map<label>, faceRegion, iter)
197  {
198  if (iter() == freeRegion)
199  {
200  iter() = useRegion;
201  }
202  }
203  }
204  }
205  }
206  }
207  }
208 }
209 
210 
211 bool Foam::combineFaces::faceNeighboursValid
212 (
213  const label cellI,
214  const Map<label>& faceRegion
215 ) const
216 {
217  if (faceRegion.size() <= 1)
218  {
219  return true;
220  }
221 
222  const cell& cFaces = mesh_.cells()[cellI];
223 
224  DynamicList<label> storage;
225 
226  // Test for face collapsing to edge since too many neighbours merged.
227  forAll(cFaces, cFaceI)
228  {
229  label faceI = cFaces[cFaceI];
230 
231  if (!faceRegion.found(faceI))
232  {
233  const labelList& fEdges = mesh_.faceEdges(faceI, storage);
234 
235  // Count number of remaining faces neighbouring faceI. This has
236  // to be 3 or more.
237 
238  // Unregioned neighbouring faces
239  DynamicList<label> neighbourFaces(cFaces.size());
240  // Regioned neighbouring faces
241  labelHashSet neighbourRegions;
242 
243  forAll(fEdges, i)
244  {
245  label edgeI = fEdges[i];
246  label nbrI = meshTools::otherFace(mesh_, cellI, faceI, edgeI);
247 
248  Map<label>::const_iterator iter = faceRegion.find(nbrI);
249 
250  if (iter == faceRegion.end())
251  {
252  if (findIndex(neighbourFaces, nbrI) == -1)
253  {
254  neighbourFaces.append(nbrI);
255  }
256  }
257  else
258  {
259  neighbourRegions.insert(iter());
260  }
261  }
262 
263  if ((neighbourFaces.size()+neighbourRegions.size()) < 3)
264  {
265  return false;
266  }
267  }
268  }
269  return true;
270 }
271 
272 
273 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
274 
275 // Construct from mesh
276 Foam::combineFaces::combineFaces
277 (
278  const polyMesh& mesh,
279  const bool undoable
280 )
281 :
282  mesh_(mesh),
283  undoable_(undoable),
284  masterFace_(0),
285  faceSetsVertices_(0),
286  savedPointLabels_(0),
287  savedPoints_(0)
288 {}
289 
290 
291 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
292 
294 (
295  const scalar featureCos,
296  const scalar minConcaveCos,
297  const labelHashSet& boundaryCells
298 ) const
299 {
300  // Lists of faces that can be merged.
301  DynamicList<labelList> allFaceSets(boundaryCells.size() / 10);
302  // Storage for on-the-fly cell-edge addressing.
303  DynamicList<label> storage;
304 
305  // On all cells regionise the faces
306  forAllConstIter(labelHashSet, boundaryCells, iter)
307  {
308  label cellI = iter.key();
309 
310  const cell& cFaces = mesh_.cells()[cellI];
311 
312  const labelList& cEdges = mesh_.cellEdges(cellI, storage);
313 
314  // Region per face
315  Map<label> faceRegion(cFaces.size());
316  regioniseFaces(featureCos, cellI, cEdges, faceRegion);
317 
318  // Now we have in faceRegion for every face the region with planar
319  // face sharing the same region. We now check whether the resulting
320  // sets cause a face
321  // - to become a set of edges since too many faces are merged.
322  // - to become convex
323 
324  if (faceNeighboursValid(cellI, faceRegion))
325  {
326  // Create region-to-faces addressing
327  Map<labelList> regionToFaces(faceRegion.size());
328 
329  forAllConstIter(Map<label>, faceRegion, iter)
330  {
331  label faceI = iter.key();
332  label region = iter();
333 
334  Map<labelList>::iterator regionFnd = regionToFaces.find(region);
335 
336  if (regionFnd != regionToFaces.end())
337  {
338  labelList& setFaces = regionFnd();
339  label sz = setFaces.size();
340  setFaces.setSize(sz+1);
341  setFaces[sz] = faceI;
342  }
343  else
344  {
345  regionToFaces.insert(region, labelList(1, faceI));
346  }
347  }
348 
349  // For every set check if it forms a valid convex face
350 
351  forAllConstIter(Map<labelList>, regionToFaces, iter)
352  {
353  // Make face out of setFaces
354  indirectPrimitivePatch bigFace
355  (
357  (
358  mesh_.faces(),
359  iter()
360  ),
361  mesh_.points()
362  );
363 
364  // Only store if -only one outside loop -which forms convex face
365  if (validFace(minConcaveCos, bigFace))
366  {
367  allFaceSets.append(iter());
368  }
369  }
370  }
371  }
372 
373  return allFaceSets.shrink();
374 }
375 
376 
378 (
379  const scalar featureCos,
380  const scalar minConcaveCos
381 ) const
382 {
383  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
384 
385  // Pick up all cells on boundary
386  labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
387 
388  forAll(patches, patchI)
389  {
390  const polyPatch& patch = patches[patchI];
391 
392  if (!patch.coupled())
393  {
394  forAll(patch, i)
395  {
396  boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
397  }
398  }
399  }
400 
401  return getMergeSets(featureCos, minConcaveCos, boundaryCells);
402 }
403 
404 
405 // Gets outside edgeloop as a face
406 // - in same order as faces
407 // - in mesh vertex labels
409 (
410  const indirectPrimitivePatch& fp
411 )
412 {
413  if (fp.edgeLoops().size() != 1)
414  {
416  (
417  "combineFaces::getOutsideFace(const indirectPrimitivePatch&)"
418  ) << "Multiple outside loops:" << fp.edgeLoops()
419  << abort(FatalError);
420  }
421 
422  // Get first boundary edge. Since guaranteed one edgeLoop when in here this
423  // edge must be on it.
424  label bEdgeI = fp.nInternalEdges();
425 
426  const edge& e = fp.edges()[bEdgeI];
427 
428  const labelList& eFaces = fp.edgeFaces()[bEdgeI];
429 
430  if (eFaces.size() != 1)
431  {
433  (
434  "combineFaces::getOutsideFace(const indirectPrimitivePatch&)"
435  ) << "boundary edge:" << bEdgeI
436  << " points:" << fp.meshPoints()[e[0]]
437  << ' ' << fp.meshPoints()[e[1]]
438  << " on indirectPrimitivePatch has " << eFaces.size()
439  << " faces using it" << abort(FatalError);
440  }
441 
442 
443  // Outside loop
444  const labelList& outsideLoop = fp.edgeLoops()[0];
445 
446 
447  // Get order of edge e in outside loop
448  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
449 
450  // edgeLoopConsistent : true = edge in same order as outsideloop
451  // false = opposite order
452  bool edgeLoopConsistent = false;
453 
454  {
455  label index0 = findIndex(outsideLoop, e[0]);
456  label index1 = findIndex(outsideLoop, e[1]);
457 
458  if (index0 == -1 || index1 == -1)
459  {
461  (
462  "combineFaces::getOutsideFace"
463  "(const indirectPrimitivePatch&)"
464  ) << "Cannot find boundary edge:" << e
465  << " points:" << fp.meshPoints()[e[0]]
466  << ' ' << fp.meshPoints()[e[1]]
467  << " in edgeLoop:" << outsideLoop << abort(FatalError);
468  }
469  else if (index1 == outsideLoop.fcIndex(index0))
470  {
471  edgeLoopConsistent = true;
472  }
473  else if (index0 == outsideLoop.fcIndex(index1))
474  {
475  edgeLoopConsistent = false;
476  }
477  else
478  {
480  (
481  "combineFaces::getOutsideFace"
482  "(const indirectPrimitivePatch&)"
483  ) << "Cannot find boundary edge:" << e
484  << " points:" << fp.meshPoints()[e[0]]
485  << ' ' << fp.meshPoints()[e[1]]
486  << " on consecutive points in edgeLoop:"
487  << outsideLoop << abort(FatalError);
488  }
489  }
490 
491 
492  // Get face in local vertices
493  const face& localF = fp.localFaces()[eFaces[0]];
494 
495  // Get order of edge in localF
496  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
497 
498  // faceEdgeConsistent : true = edge in same order as localF
499  // false = opposite order
500  bool faceEdgeConsistent = false;
501 
502  {
503  // Find edge in face.
504  label index = findIndex(fp.faceEdges()[eFaces[0]], bEdgeI);
505 
506  if (index == -1)
507  {
509  (
510  "combineFaces::getOutsideFace"
511  "(const indirectPrimitivePatch&)"
512  ) << "Cannot find boundary edge:" << e
513  << " points:" << fp.meshPoints()[e[0]]
514  << ' ' << fp.meshPoints()[e[1]]
515  << " in face:" << eFaces[0]
516  << " edges:" << fp.faceEdges()[eFaces[0]]
517  << abort(FatalError);
518  }
519  else if (localF[index] == e[0] && localF.nextLabel(index) == e[1])
520  {
521  faceEdgeConsistent = true;
522  }
523  else if (localF[index] == e[1] && localF.nextLabel(index) == e[0])
524  {
525  faceEdgeConsistent = false;
526  }
527  else
528  {
530  (
531  "combineFaces::getOutsideFace"
532  "(const indirectPrimitivePatch&)"
533  ) << "Cannot find boundary edge:" << e
534  << " points:" << fp.meshPoints()[e[0]]
535  << ' ' << fp.meshPoints()[e[1]]
536  << " in face:" << eFaces[0] << " verts:" << localF
537  << abort(FatalError);
538  }
539  }
540 
541  // Get face in mesh points.
542  face meshFace(renumber(fp.meshPoints(), outsideLoop));
543 
544  if (faceEdgeConsistent != edgeLoopConsistent)
545  {
546  reverse(meshFace);
547  }
548  return meshFace;
549 }
550 
551 
553 (
554  const labelListList& faceSets,
555  polyTopoChange& meshMod
556 )
557 {
558  if (undoable_)
559  {
560  masterFace_.setSize(faceSets.size());
561  faceSetsVertices_.setSize(faceSets.size());
562  forAll(faceSets, setI)
563  {
564  const labelList& setFaces = faceSets[setI];
565 
566  masterFace_[setI] = setFaces[0];
567  faceSetsVertices_[setI] = UIndirectList<face>
568  (
569  mesh_.faces(),
570  setFaces
571  );
572  }
573  }
574 
575  // Running count of number of faces using a point
576  labelList nPointFaces(mesh_.nPoints(), 0);
577 
578  const labelListList& pointFaces = mesh_.pointFaces();
579 
580  forAll(pointFaces, pointI)
581  {
582  nPointFaces[pointI] = pointFaces[pointI].size();
583  }
584 
585  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
586 
587  forAll(faceSets, setI)
588  {
589  const labelList& setFaces = faceSets[setI];
590 
591  // Check
592  if (debug)
593  {
594  forAll(setFaces, i)
595  {
596  label patchI = patches.whichPatch(setFaces[i]);
597 
598  if (patchI == -1 || patches[patchI].coupled())
599  {
601  (
602  "combineFaces::setRefinement"
603  "(const bool, const labelListList&"
604  ", polyTopoChange&)"
605  ) << "Can only merge non-coupled boundary faces"
606  << " but found internal or coupled face:"
607  << setFaces[i] << " in set " << setI
608  << abort(FatalError);
609  }
610  }
611  }
612 
613  // Make face out of setFaces
614  indirectPrimitivePatch bigFace
615  (
617  (
618  mesh_.faces(),
619  setFaces
620  ),
621  mesh_.points()
622  );
623 
624  // Get outside vertices (in local vertex numbering)
625  const labelListList& edgeLoops = bigFace.edgeLoops();
626 
627  if (edgeLoops.size() != 1)
628  {
630  (
631  "combineFaces::setRefinement"
632  "(const bool, const labelListList&, polyTopoChange&)"
633  ) << "Faces to-be-merged " << setFaces
634  << " do not form a single big face." << nl
635  << abort(FatalError);
636  }
637 
638 
639  // Delete all faces except master, whose face gets modified.
640 
641  // Modify master face
642  // ~~~~~~~~~~~~~~~~~~
643 
644  label masterFaceI = setFaces[0];
645 
646  // Get outside face in mesh vertex labels
647  face outsideFace(getOutsideFace(bigFace));
648 
649  label zoneID = mesh_.faceZones().whichZone(masterFaceI);
650 
651  bool zoneFlip = false;
652 
653  if (zoneID >= 0)
654  {
655  const faceZone& fZone = mesh_.faceZones()[zoneID];
656 
657  zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
658  }
659 
660  label patchI = mesh_.boundaryMesh().whichPatch(masterFaceI);
661 
662  meshMod.setAction
663  (
665  (
666  outsideFace, // modified face
667  masterFaceI, // label of face being modified
668  mesh_.faceOwner()[masterFaceI], // owner
669  -1, // neighbour
670  false, // face flip
671  patchI, // patch for face
672  false, // remove from zone
673  zoneID, // zone for face
674  zoneFlip // face flip in zone
675  )
676  );
677 
678 
679  // Delete all non-master faces
680  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
681 
682  for (label i = 1; i < setFaces.size(); i++)
683  {
684  meshMod.setAction(polyRemoveFace(setFaces[i]));
685  }
686 
687 
688  // Mark unused points for deletion
689  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
690 
691  // Remove count of all faces combined
692  forAll(setFaces, i)
693  {
694  const face& f = mesh_.faces()[setFaces[i]];
695 
696  forAll(f, fp)
697  {
698  nPointFaces[f[fp]]--;
699  }
700  }
701  // But keep count on new outside face
702  forAll(outsideFace, fp)
703  {
704  nPointFaces[outsideFace[fp]]++;
705  }
706  }
707 
708 
709  // If one or more processors use the point it needs to be kept.
711  (
712  mesh_,
713  nPointFaces,
714  plusEqOp<label>(),
715  0, // null value
716  false // no separation
717  );
718 
719  // Remove all unused points. Store position if undoable.
720  if (!undoable_)
721  {
722  forAll(nPointFaces, pointI)
723  {
724  if (nPointFaces[pointI] == 0)
725  {
726  meshMod.setAction(polyRemovePoint(pointI));
727  }
728  }
729  }
730  else
731  {
732  // Count removed points
733  label n = 0;
734  forAll(nPointFaces, pointI)
735  {
736  if (nPointFaces[pointI] == 0)
737  {
738  n++;
739  }
740  }
741 
742  savedPointLabels_.setSize(n);
743  savedPoints_.setSize(n);
744  Map<label> meshToSaved(2*n);
745 
746  // Remove points and store position
747  n = 0;
748  forAll(nPointFaces, pointI)
749  {
750  if (nPointFaces[pointI] == 0)
751  {
752  meshMod.setAction(polyRemovePoint(pointI));
753 
754  savedPointLabels_[n] = pointI;
755  savedPoints_[n] = mesh_.points()[pointI];
756 
757  meshToSaved.insert(pointI, n);
758  n++;
759  }
760  }
761 
762  // Update stored vertex labels. Negative indices index into local points
763  forAll(faceSetsVertices_, setI)
764  {
765  faceList& setFaces = faceSetsVertices_[setI];
766 
767  forAll(setFaces, i)
768  {
769  face& f = setFaces[i];
770 
771  forAll(f, fp)
772  {
773  label pointI = f[fp];
774 
775  if (nPointFaces[pointI] == 0)
776  {
777  f[fp] = -meshToSaved[pointI]-1;
778  }
779  }
780  }
781  }
782  }
783 }
784 
785 
787 {
788  if (undoable_)
789  {
790  // Master face just renumbering of point labels
791  inplaceRenumber(map.reverseFaceMap(), masterFace_);
792 
793  // Stored faces refer to backed-up vertices (not changed)
794  // and normal vertices (need to be renumbered)
795  forAll(faceSetsVertices_, setI)
796  {
797  faceList& faces = faceSetsVertices_[setI];
798 
799  forAll(faces, i)
800  {
801  // Note: faces[i] can have negative elements.
802  face& f = faces[i];
803 
804  forAll(f, fp)
805  {
806  label pointI = f[fp];
807 
808  if (pointI >= 0)
809  {
810  f[fp] = map.reversePointMap()[pointI];
811 
812  if (f[fp] < 0)
813  {
815  (
816  "combineFaces::updateMesh"
817  "(const mapPolyMesh&)"
818  ) << "In set " << setI << " at position " << i
819  << " with master face "
820  << masterFace_[setI] << nl
821  << "the points of the slave face " << faces[i]
822  << " don't exist anymore."
823  << abort(FatalError);
824  }
825  }
826  }
827  }
828  }
829  }
830 }
831 
832 
833 // Note that faces on coupled patches are never combined (or at least
834 // getMergeSets never picks them up. Thus the points on them will never get
835 // deleted since still used by the coupled faces. So never a risk of undoing
836 // things (faces/points) on coupled patches.
838 (
839  const labelList& masterFaces,
840  polyTopoChange& meshMod,
841  Map<label>& restoredPoints,
842  Map<label>& restoredFaces,
843  Map<label>& restoredCells
844 )
845 {
846  if (!undoable_)
847  {
849  (
850  "combineFaces::setUnrefinement"
851  "(const labelList&, polyTopoChange&"
852  ", Map<label>&, Map<label>&, Map<label>&)"
853  ) << "Can only call setUnrefinement if constructed with"
854  << " unrefinement capability." << exit(FatalError);
855  }
856 
857 
858  // Restored points
859  labelList addedPoints(savedPoints_.size(), -1);
860 
861  // Invert set-to-master-face
862  Map<label> masterToSet(masterFace_.size());
863 
864  forAll(masterFace_, setI)
865  {
866  if (masterFace_[setI] >= 0)
867  {
868  masterToSet.insert(masterFace_[setI], setI);
869  }
870  }
871 
872  forAll(masterFaces, i)
873  {
874  label masterFaceI = masterFaces[i];
875 
876  Map<label>::const_iterator iter = masterToSet.find(masterFaceI);
877 
878  if (iter == masterToSet.end())
879  {
881  (
882  "combineFaces::setUnrefinement"
883  "(const labelList&, polyTopoChange&"
884  ", Map<label>&, Map<label>&, Map<label>&)"
885  ) << "Master face " << masterFaceI
886  << " is not the master of one of the merge sets"
887  << " or has already been merged"
888  << abort(FatalError);
889  }
890 
891  label setI = iter();
892 
893 
894  // Update faces of the merge setI for reintroduced vertices
895  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
896 
897  faceList& faces = faceSetsVertices_[setI];
898 
899  if (faces.empty())
900  {
902  (
903  "combineFaces::setUnrefinement"
904  "(const labelList&, polyTopoChange&"
905  ", Map<label>&, Map<label>&, Map<label>&)"
906  ) << "Set " << setI << " with master face " << masterFaceI
907  << " has already been merged." << abort(FatalError);
908  }
909 
910  forAll(faces, i)
911  {
912  face& f = faces[i];
913 
914  forAll(f, fp)
915  {
916  label pointI = f[fp];
917 
918  if (pointI < 0)
919  {
920  label localI = -pointI-1;
921 
922  if (addedPoints[localI] == -1)
923  {
924  // First occurrence of saved point. Reintroduce point
925  addedPoints[localI] = meshMod.setAction
926  (
928  (
929  savedPoints_[localI], // point
930  -1, // master point
931  -1, // zone for point
932  true // supports a cell
933  )
934  );
935  restoredPoints.insert
936  (
937  addedPoints[localI], // current point label
938  savedPointLabels_[localI] // point label when it
939  // was stored
940  );
941  }
942  f[fp] = addedPoints[localI];
943  }
944  }
945  }
946 
947 
948  // Restore
949  // ~~~~~~~
950 
951  label own = mesh_.faceOwner()[masterFaceI];
952  label zoneID = mesh_.faceZones().whichZone(masterFaceI);
953  bool zoneFlip = false;
954  if (zoneID >= 0)
955  {
956  const faceZone& fZone = mesh_.faceZones()[zoneID];
957  zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
958  }
959  label patchI = mesh_.boundaryMesh().whichPatch(masterFaceI);
960 
961  if (mesh_.boundaryMesh()[patchI].coupled())
962  {
964  (
965  "combineFaces::setUnrefinement"
966  "(const labelList&, polyTopoChange&"
967  ", Map<label>&, Map<label>&, Map<label>&)"
968  ) << "Master face " << masterFaceI << " is on coupled patch "
969  << mesh_.boundaryMesh()[patchI].name()
970  << abort(FatalError);
971  }
972 
973  //Pout<< "Restoring new master face " << masterFaceI
974  // << " to vertices " << faces[0] << endl;
975 
976  // Modify the master face.
977  meshMod.setAction
978  (
980  (
981  faces[0], // original face
982  masterFaceI, // label of face
983  own, // owner
984  -1, // neighbour
985  false, // face flip
986  patchI, // patch for face
987  false, // remove from zone
988  zoneID, // zone for face
989  zoneFlip // face flip in zone
990  )
991  );
992 
993  // Add the previously removed faces
994  for (label i = 1; i < faces.size(); i++)
995  {
996  //Pout<< "Restoring removed face with vertices " << faces[i]
997  // << endl;
998 
999  meshMod.setAction
1000  (
1001  polyAddFace
1002  (
1003  faces[i], // vertices
1004  own, // owner,
1005  -1, // neighbour,
1006  -1, // masterPointID,
1007  -1, // masterEdgeID,
1008  masterFaceI, // masterFaceID,
1009  false, // flipFaceFlux,
1010  patchI, // patchID,
1011  zoneID, // zoneID,
1012  zoneFlip // zoneFlip
1013  )
1014  );
1015  }
1016 
1017  // Clear out restored set
1018  faceSetsVertices_[setI].clear();
1019  masterFace_[setI] = -1;
1020  }
1021 }
1022 
1023 
1024 // ************************ vim: set sw=4 sts=4 et: ************************ //