FreeFOAM The Cross-Platform CFD Toolkit
removeFaces.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 "removeFaces.H"
27 #include <OpenFOAM/polyMesh.H>
28 #include "polyTopoChange.H"
29 #include <meshTools/meshTools.H>
34 #include <OpenFOAM/syncTools.H>
35 #include <OpenFOAM/OFstream.H>
37 #include <OpenFOAM/Time.H>
38 #include <meshTools/faceSet.H>
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
42 namespace Foam
43 {
44 
45 defineTypeNameAndDebug(removeFaces, 0);
46 
47 }
48 
49 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
50 
51 // Changes region of connected set of cells. Can be recursive since hopefully
52 // only small area of faces removed in one go.
53 void Foam::removeFaces::changeCellRegion
54 (
55  const label cellI,
56  const label oldRegion,
57  const label newRegion,
58  labelList& cellRegion
59 ) const
60 {
61  if (cellRegion[cellI] == oldRegion)
62  {
63  cellRegion[cellI] = newRegion;
64 
65  // Step to neighbouring cells
66 
67  const labelList& cCells = mesh_.cellCells()[cellI];
68 
69  forAll(cCells, i)
70  {
71  changeCellRegion(cCells[i], oldRegion, newRegion, cellRegion);
72  }
73  }
74 }
75 
76 
77 // Changes region of connected set of faces. Returns number of changed faces.
78 Foam::label Foam::removeFaces::changeFaceRegion
79 (
80  const labelList& cellRegion,
81  const boolList& removedFace,
82  const labelList& nFacesPerEdge,
83  const label faceI,
84  const label newRegion,
85  const labelList& fEdges,
86  labelList& faceRegion
87 ) const
88 {
89  label nChanged = 0;
90 
91  if (faceRegion[faceI] == -1 && !removedFace[faceI])
92  {
93  faceRegion[faceI] = newRegion;
94 
95  nChanged = 1;
96 
97  // Storage for on-the-fly addressing
98  DynamicList<label> fe;
99  DynamicList<label> ef;
100 
101  // Step to neighbouring faces across edges that will get removed
102  forAll(fEdges, i)
103  {
104  label edgeI = fEdges[i];
105 
106  if (nFacesPerEdge[edgeI] >= 0 && nFacesPerEdge[edgeI] <= 2)
107  {
108  const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
109 
110  forAll(eFaces, j)
111  {
112  label nbrFaceI = eFaces[j];
113 
114  const labelList& fEdges1 = mesh_.faceEdges(nbrFaceI, fe);
115 
116  nChanged += changeFaceRegion
117  (
118  cellRegion,
119  removedFace,
120  nFacesPerEdge,
121  nbrFaceI,
122  newRegion,
123  fEdges1,
124  faceRegion
125  );
126  }
127  }
128  }
129  }
130  return nChanged;
131 }
132 
133 
134 // Mark all faces affected in any way by
135 // - removal of cells
136 // - removal of faces
137 // - removal of edges
138 // - removal of points
139 Foam::boolList Foam::removeFaces::getFacesAffected
140 (
141  const labelList& cellRegion,
142  const labelList& cellRegionMaster,
143  const labelList& facesToRemove,
144  const labelHashSet& edgesToRemove,
145  const labelHashSet& pointsToRemove
146 ) const
147 {
148  boolList affectedFace(mesh_.nFaces(), false);
149 
150  // Mark faces affected by removal of cells
151  forAll(cellRegion, cellI)
152  {
153  label region = cellRegion[cellI];
154 
155  if (region != -1 && (cellI != cellRegionMaster[region]))
156  {
157  const labelList& cFaces = mesh_.cells()[cellI];
158 
159  forAll(cFaces, cFaceI)
160  {
161  affectedFace[cFaces[cFaceI]] = true;
162  }
163  }
164  }
165 
166  // Mark faces affected by removal of face.
167  forAll(facesToRemove, i)
168  {
169  affectedFace[facesToRemove[i]] = true;
170  }
171 
172  // Mark faces affected by removal of edges
173  forAllConstIter(labelHashSet, edgesToRemove, iter)
174  {
175  const labelList& eFaces = mesh_.edgeFaces(iter.key());
176 
177  forAll(eFaces, eFaceI)
178  {
179  affectedFace[eFaces[eFaceI]] = true;
180  }
181  }
182 
183  // Mark faces affected by removal of points
184  forAllConstIter(labelHashSet, pointsToRemove, iter)
185  {
186  label pointI = iter.key();
187 
188  const labelList& pFaces = mesh_.pointFaces()[pointI];
189 
190  forAll(pFaces, pFaceI)
191  {
192  affectedFace[pFaces[pFaceI]] = true;
193  }
194  }
195  return affectedFace;
196 }
197 
198 
200 (
201  const indirectPrimitivePatch& fp,
202  const fileName& fName
203 )
204 {
205  OFstream str(fName);
206  Pout<< "removeFaces::writeOBJ : Writing faces to file "
207  << str.name() << endl;
208 
209  const pointField& localPoints = fp.localPoints();
210 
211  forAll(localPoints, i)
212  {
213  meshTools::writeOBJ(str, localPoints[i]);
214  }
215 
216  const faceList& localFaces = fp.localFaces();
217 
218  forAll(localFaces, i)
219  {
220  const face& f = localFaces[i];
221 
222  str<< 'f';
223 
224  forAll(f, fp)
225  {
226  str<< ' ' << f[fp]+1;
227  }
228  str<< nl;
229  }
230 }
231 
232 
233 // Inserts commands to merge faceLabels into one face.
234 void Foam::removeFaces::mergeFaces
235 (
236  const labelList& cellRegion,
237  const labelList& cellRegionMaster,
238  const labelHashSet& pointsToRemove,
239  const labelList& faceLabels,
240  polyTopoChange& meshMod
241 ) const
242 {
243  // Construct addressing engine from faceLabels (in order of faceLabels as
244  // well)
246  (
247  IndirectList<face>
248  (
249  mesh_.faces(),
250  faceLabels
251  ),
252  mesh_.points()
253  );
254 
255  // Get outside vertices (in local vertex numbering)
256 
257  if (fp.edgeLoops().size() != 1)
258  {
259  writeOBJ(fp, mesh_.time().path()/"facesToBeMerged.obj");
260  FatalErrorIn("removeFaces::mergeFaces")
261  << "Cannot merge faces " << faceLabels
262  << " into single face since outside vertices " << fp.edgeLoops()
263  << " do not form single loop but form " << fp.edgeLoops().size()
264  << " loops instead." << abort(FatalError);
265  }
266 
267  const labelList& edgeLoop = fp.edgeLoops()[0];
268 
269  // Get outside vertices in order of one of the faces in faceLabels.
270  // (this becomes the master face)
271  // Find the first face that uses edgeLoop[0] and edgeLoop[1] as consecutive
272  // vertices.
273 
274  label masterIndex = -1;
275  bool reverseLoop = false;
276 
277  const labelList& pFaces = fp.pointFaces()[edgeLoop[0]];
278 
279  // Find face among pFaces which uses edgeLoop[1]
280  forAll(pFaces, i)
281  {
282  label faceI = pFaces[i];
283 
284  const face& f = fp.localFaces()[faceI];
285 
286  label index1 = findIndex(f, edgeLoop[1]);
287 
288  if (index1 != -1)
289  {
290  // Check whether consecutive to edgeLoop[0]
291  label index0 = findIndex(f, edgeLoop[0]);
292 
293  if (index0 != -1)
294  {
295  if (index1 == f.fcIndex(index0))
296  {
297  masterIndex = faceI;
298  reverseLoop = false;
299  break;
300  }
301  else if (index1 == f.rcIndex(index0))
302  {
303  masterIndex = faceI;
304  reverseLoop = true;
305  break;
306  }
307  }
308  }
309  }
310 
311  if (masterIndex == -1)
312  {
313  writeOBJ(fp, mesh_.time().path()/"facesToBeMerged.obj");
314  FatalErrorIn("removeFaces::mergeFaces")
315  << "Problem" << abort(FatalError);
316  }
317 
318 
319  // Modify the master face.
320  // ~~~~~~~~~~~~~~~~~~~~~~~
321 
322  // Modify first face.
323  label faceI = faceLabels[masterIndex];
324 
325  label own = mesh_.faceOwner()[faceI];
326 
327  if (cellRegion[own] != -1)
328  {
329  own = cellRegionMaster[cellRegion[own]];
330  }
331 
332  label patchID, zoneID, zoneFlip;
333 
334  getFaceInfo(faceI, patchID, zoneID, zoneFlip);
335 
336  label nei = -1;
337 
338  if (mesh_.isInternalFace(faceI))
339  {
340  nei = mesh_.faceNeighbour()[faceI];
341 
342  if (cellRegion[nei] != -1)
343  {
344  nei = cellRegionMaster[cellRegion[nei]];
345  }
346  }
347 
348 
349  DynamicList<label> faceVerts(edgeLoop.size());
350 
351  forAll(edgeLoop, i)
352  {
353  label pointI = fp.meshPoints()[edgeLoop[i]];
354 
355  if (pointsToRemove.found(pointI))
356  {
357  //Pout<< "**Removing point " << pointI << " from "
358  // << edgeLoop << endl;
359  }
360  else
361  {
362  faceVerts.append(pointI);
363  }
364  }
365 
366  face mergedFace;
367  mergedFace.transfer(faceVerts);
368 
369  if (reverseLoop)
370  {
371  reverse(mergedFace);
372  }
373 
374  //{
375  // Pout<< "Modifying masterface " << faceI
376  // << " from faces:" << faceLabels
377  // << " old verts:" << UIndirectList<face>(mesh_.faces(), faceLabels)
378  // << " for new verts:"
379  // << mergedFace
380  // << " possibly new owner " << own
381  // << " or new nei " << nei
382  // << endl;
383  //}
384 
385  modFace
386  (
387  mergedFace, // modified face
388  faceI, // label of face being modified
389  own, // owner
390  nei, // neighbour
391  false, // face flip
392  patchID, // patch for face
393  false, // remove from zone
394  zoneID, // zone for face
395  zoneFlip, // face flip in zone
396 
397  meshMod
398  );
399 
400 
401  // Remove all but master face.
402  forAll(faceLabels, patchFaceI)
403  {
404  if (patchFaceI != masterIndex)
405  {
406  //Pout<< "Removing face " << faceLabels[patchFaceI] << endl;
407 
408  meshMod.setAction(polyRemoveFace(faceLabels[patchFaceI], faceI));
409  }
410  }
411 }
412 
413 
414 // Get patch, zone info for faceI
415 void Foam::removeFaces::getFaceInfo
416 (
417  const label faceI,
418 
419  label& patchID,
420  label& zoneID,
421  label& zoneFlip
422 ) const
423 {
424  patchID = -1;
425 
426  if (!mesh_.isInternalFace(faceI))
427  {
428  patchID = mesh_.boundaryMesh().whichPatch(faceI);
429  }
430 
431  zoneID = mesh_.faceZones().whichZone(faceI);
432 
433  zoneFlip = false;
434 
435  if (zoneID >= 0)
436  {
437  const faceZone& fZone = mesh_.faceZones()[zoneID];
438 
439  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
440  }
441 }
442 
443 
444 // Return face with all pointsToRemove removed.
445 Foam::face Foam::removeFaces::filterFace
446 (
447  const labelHashSet& pointsToRemove,
448  const label faceI
449 ) const
450 {
451  const face& f = mesh_.faces()[faceI];
452 
453  labelList newFace(f.size(), -1);
454 
455  label newFp = 0;
456 
457  forAll(f, fp)
458  {
459  label vertI = f[fp];
460 
461  if (!pointsToRemove.found(vertI))
462  {
463  newFace[newFp++] = vertI;
464  }
465  }
466 
467  newFace.setSize(newFp);
468 
469  return face(newFace);
470 }
471 
472 
473 // Wrapper for meshMod.modifyFace. Reverses face if own>nei.
474 void Foam::removeFaces::modFace
475 (
476  const face& f,
477  const label masterFaceID,
478  const label own,
479  const label nei,
480  const bool flipFaceFlux,
481  const label newPatchID,
482  const bool removeFromZone,
483  const label zoneID,
484  const bool zoneFlip,
485 
486  polyTopoChange& meshMod
487 ) const
488 {
489  if ((nei == -1) || (own < nei))
490  {
491 // if (debug)
492 // {
493 // Pout<< "ModifyFace (unreversed) :"
494 // << " faceI:" << masterFaceID
495 // << " f:" << f
496 // << " own:" << own
497 // << " nei:" << nei
498 // << " flipFaceFlux:" << flipFaceFlux
499 // << " newPatchID:" << newPatchID
500 // << " removeFromZone:" << removeFromZone
501 // << " zoneID:" << zoneID
502 // << " zoneFlip:" << zoneFlip
503 // << endl;
504 // }
505 
506  meshMod.setAction
507  (
508  polyModifyFace
509  (
510  f, // modified face
511  masterFaceID, // label of face being modified
512  own, // owner
513  nei, // neighbour
514  flipFaceFlux, // face flip
515  newPatchID, // patch for face
516  removeFromZone, // remove from zone
517  zoneID, // zone for face
518  zoneFlip // face flip in zone
519  )
520  );
521  }
522  else
523  {
524 // if (debug)
525 // {
526 // Pout<< "ModifyFace (!reversed) :"
527 // << " faceI:" << masterFaceID
528 // << " f:" << f.reverseFace()
529 // << " own:" << nei
530 // << " nei:" << own
531 // << " flipFaceFlux:" << flipFaceFlux
532 // << " newPatchID:" << newPatchID
533 // << " removeFromZone:" << removeFromZone
534 // << " zoneID:" << zoneID
535 // << " zoneFlip:" << zoneFlip
536 // << endl;
537 // }
538 
539  meshMod.setAction
540  (
541  polyModifyFace
542  (
543  f.reverseFace(),// modified face
544  masterFaceID, // label of face being modified
545  nei, // owner
546  own, // neighbour
547  flipFaceFlux, // face flip
548  newPatchID, // patch for face
549  removeFromZone, // remove from zone
550  zoneID, // zone for face
551  zoneFlip // face flip in zone
552  )
553  );
554  }
555 }
556 
557 
558 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
559 
560 // Construct from mesh
561 Foam::removeFaces::removeFaces
562 (
563  const polyMesh& mesh,
564  const scalar minCos
565 )
566 :
567  mesh_(mesh),
568  minCos_(minCos)
569 {}
570 
571 
572 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
573 
574 // Removing face connects cells. This function works out a consistent set of
575 // cell regions.
576 // - returns faces to remove. Can be extended with additional faces
577 // (if owner would become neighbour)
578 // - sets cellRegion to -1 or to region number
579 // - regionMaster contains for every region number a master cell.
581 (
582  const labelList& facesToRemove,
583  labelList& cellRegion,
584  labelList& regionMaster,
585  labelList& newFacesToRemove
586 ) const
587 {
588  const labelList& faceOwner = mesh_.faceOwner();
589  const labelList& faceNeighbour = mesh_.faceNeighbour();
590 
591  cellRegion.setSize(mesh_.nCells());
592  cellRegion = -1;
593 
594  regionMaster.setSize(mesh_.nCells());
595  regionMaster = -1;
596 
597  label nRegions = 0;
598 
599  forAll(facesToRemove, i)
600  {
601  label faceI = facesToRemove[i];
602 
603  if (!mesh_.isInternalFace(faceI))
604  {
606  (
607  "removeFaces::compatibleRemoves(const labelList&"
608  ", labelList&, labelList&, labelList&)"
609  ) << "Not internal face:" << faceI << abort(FatalError);
610  }
611 
612 
613  label own = faceOwner[faceI];
614  label nei = faceNeighbour[faceI];
615 
616  label region0 = cellRegion[own];
617  label region1 = cellRegion[nei];
618 
619  if (region0 == -1)
620  {
621  if (region1 == -1)
622  {
623  // Create new region
624  cellRegion[own] = nRegions;
625  cellRegion[nei] = nRegions;
626 
627  // Make owner (lowest numbered!) the master of the region
628  regionMaster[nRegions] = own;
629  nRegions++;
630  }
631  else
632  {
633  // Add owner to neighbour region
634  cellRegion[own] = region1;
635  // See if owner becomes the master of the region
636  regionMaster[region1] = min(own, regionMaster[region1]);
637  }
638  }
639  else
640  {
641  if (region1 == -1)
642  {
643  // Add neighbour to owner region
644  cellRegion[nei] = region0;
645  // nei is higher numbered than own so guaranteed not lower
646  // than master of region0.
647  }
648  else if (region0 != region1)
649  {
650  // Both have regions. Keep lowest numbered region and master.
651  label freedRegion = -1;
652  label keptRegion = -1;
653 
654  if (region0 < region1)
655  {
656  changeCellRegion
657  (
658  nei,
659  region1, // old region
660  region0, // new region
661  cellRegion
662  );
663 
664  keptRegion = region0;
665  freedRegion = region1;
666  }
667  else if (region1 < region0)
668  {
669  changeCellRegion
670  (
671  own,
672  region0, // old region
673  region1, // new region
674  cellRegion
675  );
676 
677  keptRegion = region1;
678  freedRegion = region0;
679  }
680 
681  label master0 = regionMaster[region0];
682  label master1 = regionMaster[region1];
683 
684  regionMaster[freedRegion] = -1;
685  regionMaster[keptRegion] = min(master0, master1);
686  }
687  }
688  }
689 
690  regionMaster.setSize(nRegions);
691 
692 
693  // Various checks
694  // - master is lowest numbered in any region
695  // - regions have more than 1 cell
696  {
697  labelList nCells(regionMaster.size(), 0);
698 
699  forAll(cellRegion, cellI)
700  {
701  label r = cellRegion[cellI];
702 
703  if (r != -1)
704  {
705  nCells[r]++;
706 
707  if (cellI < regionMaster[r])
708  {
710  (
711  "removeFaces::compatibleRemoves(const labelList&"
712  ", labelList&, labelList&, labelList&)"
713  ) << "Not lowest numbered : cell:" << cellI
714  << " region:" << r
715  << " regionmaster:" << regionMaster[r]
716  << abort(FatalError);
717  }
718  }
719  }
720 
721  forAll(nCells, region)
722  {
723  if (nCells[region] == 1)
724  {
726  (
727  "removeFaces::compatibleRemoves(const labelList&"
728  ", labelList&, labelList&, labelList&)"
729  ) << "Region " << region
730  << " has only " << nCells[region] << " cells in it"
731  << abort(FatalError);
732  }
733  }
734  }
735 
736 
737  // Count number of used regions
738  label nUsedRegions = 0;
739 
740  forAll(regionMaster, i)
741  {
742  if (regionMaster[i] != -1)
743  {
744  nUsedRegions++;
745  }
746  }
747 
748  // Recreate facesToRemove to be consistent with the cellRegions.
749  DynamicList<label> allFacesToRemove(facesToRemove.size());
750 
751  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
752  {
753  label own = faceOwner[faceI];
754  label nei = faceNeighbour[faceI];
755 
756  if (cellRegion[own] != -1 && cellRegion[own] == cellRegion[nei])
757  {
758  // Both will become the same cell so add face to list of faces
759  // to be removed.
760  allFacesToRemove.append(faceI);
761  }
762  }
763 
764  newFacesToRemove.transfer(allFacesToRemove);
765 
766  return nUsedRegions;
767 }
768 
769 
771 (
772  const labelList& faceLabels,
773  const labelList& cellRegion,
774  const labelList& cellRegionMaster,
775  polyTopoChange& meshMod
776 ) const
777 {
778  if (debug)
779  {
780  faceSet facesToRemove(mesh_, "facesToRemove", faceLabels);
781  Pout<< "Writing faces to remove to faceSet " << facesToRemove.name()
782  << endl;
783  facesToRemove.write();
784  }
785 
786  // Make map of all faces to be removed
787  boolList removedFace(mesh_.nFaces(), false);
788 
789  forAll(faceLabels, i)
790  {
791  label faceI = faceLabels[i];
792 
793  if (!mesh_.isInternalFace(faceI))
794  {
796  (
797  "removeFaces::setRefinement(const labelList&"
798  ", const labelList&, const labelList&, polyTopoChange&)"
799  ) << "Face to remove is not internal face:" << faceI
800  << abort(FatalError);
801  }
802 
803  removedFace[faceI] = true;
804  }
805 
806 
807  // Edges to be removed
808  // ~~~~~~~~~~~~~~~~~~~
809 
810 
811  // Edges to remove
812  labelHashSet edgesToRemove(faceLabels.size());
813 
814  // Per face the region it is in. -1 for removed faces, -2 for regions
815  // consisting of single face only.
816  labelList faceRegion(mesh_.nFaces(), -1);
817 
818  // Number of connected face regions
819  label nRegions = 0;
820 
821  // Storage for on-the-fly addressing
824 
825 
826  {
827  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
828 
829  // Usage of edges by non-removed faces.
830  // See below about initialization.
831  labelList nFacesPerEdge(mesh_.nEdges(), -1);
832 
833  // Count usage of edges by non-removed faces.
834  forAll(faceLabels, i)
835  {
836  label faceI = faceLabels[i];
837 
838  const labelList& fEdges = mesh_.faceEdges(faceI, fe);
839 
840  forAll(fEdges, i)
841  {
842  label edgeI = fEdges[i];
843 
844  if (nFacesPerEdge[edgeI] == -1)
845  {
846  nFacesPerEdge[edgeI] = mesh_.edgeFaces(edgeI, ef).size()-1;
847  }
848  else
849  {
850  nFacesPerEdge[edgeI]--;
851  }
852  }
853  }
854 
855  // Count usage for edges not on faces-to-be-removed.
856  // Note that this only needs to be done for possibly coupled edges
857  // so we could choose to loop only over boundary faces and use faceEdges
858  // of those.
859 
860  forAll(mesh_.edges(), edgeI)
861  {
862  if (nFacesPerEdge[edgeI] == -1)
863  {
864  // Edge not yet handled in loop above so is not used by any
865  // face to be removed.
866 
867  const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
868 
869  if (eFaces.size() > 2)
870  {
871  nFacesPerEdge[edgeI] = eFaces.size();
872  }
873  else if (eFaces.size() == 2)
874  {
875  // nFacesPerEdge already -1 so do nothing.
876  }
877  else
878  {
879  const edge& e = mesh_.edges()[edgeI];
880 
881  FatalErrorIn("removeFaces::setRefinement")
882  << "Problem : edge has too few face neighbours:"
883  << eFaces << endl
884  << "edge:" << edgeI
885  << " vertices:" << e
886  << " coords:" << mesh_.points()[e[0]]
887  << mesh_.points()[e[1]]
888  << abort(FatalError);
889  }
890  }
891  }
892 
893 
894 
895  if (debug)
896  {
897  OFstream str(mesh_.time().path()/"edgesWithTwoFaces.obj");
898  Pout<< "Dumping edgesWithTwoFaces to " << str.name() << endl;
899  label vertI = 0;
900 
901  forAll(nFacesPerEdge, edgeI)
902  {
903  if (nFacesPerEdge[edgeI] == 2)
904  {
905  // Edge will get removed.
906  const edge& e = mesh_.edges()[edgeI];
907 
908  meshTools::writeOBJ(str, mesh_.points()[e[0]]);
909  vertI++;
910  meshTools::writeOBJ(str, mesh_.points()[e[1]]);
911  vertI++;
912  str<< "l " << vertI-1 << ' ' << vertI << nl;
913  }
914  }
915  }
916 
917 
918  // Now all unaffected edges will have labelMax, all affected edges the
919  // number of unremoved faces.
920 
921  // Filter for edges inbetween two remaining boundary faces that
922  // make too big an angle.
923  forAll(nFacesPerEdge, edgeI)
924  {
925  if (nFacesPerEdge[edgeI] == 2)
926  {
927  // See if they are two boundary faces
928  label f0 = -1;
929  label f1 = -1;
930 
931  const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
932 
933  forAll(eFaces, i)
934  {
935  label faceI = eFaces[i];
936 
937  if (!removedFace[faceI] && !mesh_.isInternalFace(faceI))
938  {
939  if (f0 == -1)
940  {
941  f0 = faceI;
942  }
943  else
944  {
945  f1 = faceI;
946  break;
947  }
948  }
949  }
950 
951  if (f0 != -1 && f1 != -1)
952  {
953  // Edge has two boundary faces remaining.
954  // See if should be merged.
955 
956  label patch0 = patches.whichPatch(f0);
957  label patch1 = patches.whichPatch(f1);
958 
959  if (patch0 != patch1)
960  {
961  // Different patches. Do not merge edge.
962  WarningIn("removeFaces::setRefinement")
963  << "not merging faces " << f0 << " and "
964  << f1 << " across patch boundary edge " << edgeI
965  << endl;
966 
967  // Mark so it gets preserved
968  nFacesPerEdge[edgeI] = 3;
969  }
970  else if (minCos_ < 1 && minCos_ > -1)
971  {
972  const polyPatch& pp0 = patches[patch0];
973  const vectorField& n0 = pp0.faceNormals();
974 
975  if
976  (
977  mag
978  (
979  n0[f0 - pp0.start()]
980  & n0[f1 - pp0.start()]
981  )
982  < minCos_
983  )
984  {
985  WarningIn("removeFaces::setRefinement")
986  << "not merging faces " << f0 << " and "
987  << f1 << " across edge " << edgeI
988  << endl;
989 
990  // Angle between two remaining faces too large.
991  // Mark so it gets preserved
992  nFacesPerEdge[edgeI] = 3;
993  }
994  }
995  }
996  else if (f0 != -1 || f1 != -1)
997  {
998  const edge& e = mesh_.edges()[edgeI];
999 
1000  // Only found one boundary face. Problem.
1001  FatalErrorIn("removeFaces::setRefinement")
1002  << "Problem : edge would have one boundary face"
1003  << " and one internal face using it." << endl
1004  << "Your remove pattern is probably incorrect." << endl
1005  << "edge:" << edgeI
1006  << " nFaces:" << nFacesPerEdge[edgeI]
1007  << " vertices:" << e
1008  << " coords:" << mesh_.points()[e[0]]
1009  << mesh_.points()[e[1]]
1010  << " face0:" << f0
1011  << " face1:" << f1
1012  << abort(FatalError);
1013  }
1014  }
1015  }
1016 
1017 
1018 
1019  // Check locally (before synchronizing) for strangeness
1020  forAll(nFacesPerEdge, edgeI)
1021  {
1022  if (nFacesPerEdge[edgeI] == 1)
1023  {
1024  const edge& e = mesh_.edges()[edgeI];
1025 
1026  FatalErrorIn("removeFaces::setRefinement")
1027  << "Problem : edge would get 1 face using it only"
1028  << " edge:" << edgeI
1029  << " nFaces:" << nFacesPerEdge[edgeI]
1030  << " vertices:" << e
1031  << " coords:" << mesh_.points()[e[0]]
1032  << ' ' << mesh_.points()[e[1]]
1033  << abort(FatalError);
1034  }
1035  // Could check here for boundary edge with <=1 faces remaining.
1036  }
1037 
1038 
1039  // Synchronize edge usage. This is to make sure that both sides remove
1040  // (or not remove) an edge on the boundary at the same time.
1041  //
1042  // Coupled edges (edge0, edge1 are opposite each other)
1043  // a. edge not on face to be removed, edge has >= 3 faces
1044  // b. ,, edge has 2 faces
1045  // c. edge has >= 3 remaining faces
1046  // d. edge has 2 remaining faces (assume angle>minCos already handled)
1047  //
1048  // - a + a: do not remove edge
1049  // - a + b: do not remove edge
1050  // - a + c: do not remove edge
1051  // - a + d: do not remove edge
1052  //
1053  // - b + b: do not remove edge
1054  // - b + c: do not remove edge
1055  // - b + d: remove edge
1056  //
1057  // - c + c: do not remove edge
1058  // - c + d: do not remove edge
1059  // - d + d: remove edge
1060  //
1061  //
1062  // So code situation a. with >= 3
1063  // b. with -1
1064  // c. with >=3
1065  // d. with 2
1066  // then do max and check result.
1067  //
1068  // a+a : max(3,3) = 3. do not remove
1069  // a+b : max(3,-1) = 3. do not remove
1070  // a+c : max(3,3) = 3. do not remove
1071  // a+d : max(3,2) = 3. do not remove
1072  //
1073  // b+b : max(-1,-1) = -1. do not remove
1074  // b+c : max(-1,3) = 3. do not remove
1075  // b+d : max(-1,2) = 2. remove
1076  //
1077  // c+c : max(3,3) = 3. do not remove
1078  // c+d : max(3,2) = 3. do not remove
1079  //
1080  // d+d : max(2,2) = 2. remove
1081 
1083  (
1084  mesh_,
1085  nFacesPerEdge,
1086  maxEqOp<label>(),
1087  labelMin, // guaranteed to be overridden by maxEqOp
1088  false // no separation
1089  );
1090 
1091  // Convert to labelHashSet
1092  forAll(nFacesPerEdge, edgeI)
1093  {
1094  if (nFacesPerEdge[edgeI] == 0)
1095  {
1096  // 0: edge not used anymore.
1097  edgesToRemove.insert(edgeI);
1098  }
1099  else if (nFacesPerEdge[edgeI] == 1)
1100  {
1101  // 1: illegal. Tested above.
1102  }
1103  else if (nFacesPerEdge[edgeI] == 2)
1104  {
1105  // 2: merge faces.
1106  edgesToRemove.insert(edgeI);
1107  }
1108  }
1109 
1110  if (debug)
1111  {
1112  OFstream str(mesh_.time().path()/"edgesToRemove.obj");
1113  Pout<< "Dumping edgesToRemove to " << str.name() << endl;
1114  label vertI = 0;
1115 
1116  forAllConstIter(labelHashSet, edgesToRemove, iter)
1117  {
1118  // Edge will get removed.
1119  const edge& e = mesh_.edges()[iter.key()];
1120 
1121  meshTools::writeOBJ(str, mesh_.points()[e[0]]);
1122  vertI++;
1123  meshTools::writeOBJ(str, mesh_.points()[e[1]]);
1124  vertI++;
1125  str<< "l " << vertI-1 << ' ' << vertI << nl;
1126  }
1127  }
1128 
1129 
1130  // Walk to fill faceRegion with faces that will be connected across
1131  // edges that will be removed.
1132 
1133  label startFaceI = 0;
1134 
1135  while (true)
1136  {
1137  // Find unset region.
1138  for (; startFaceI < mesh_.nFaces(); startFaceI++)
1139  {
1140  if (faceRegion[startFaceI] == -1 && !removedFace[startFaceI])
1141  {
1142  break;
1143  }
1144  }
1145 
1146  if (startFaceI == mesh_.nFaces())
1147  {
1148  break;
1149  }
1150 
1151  // Start walking face-edge-face, crossing edges that will get
1152  // removed. Every thus connected region will get single region
1153  // number.
1154  label nRegion = changeFaceRegion
1155  (
1156  cellRegion,
1157  removedFace,
1158  nFacesPerEdge,
1159  startFaceI,
1160  nRegions,
1161  mesh_.faceEdges(startFaceI, fe),
1162  faceRegion
1163  );
1164 
1165  if (nRegion < 1)
1166  {
1167  FatalErrorIn("setRefinement") << "Problem" << abort(FatalError);
1168  }
1169  else if (nRegion == 1)
1170  {
1171  // Reset face to be single region.
1172  faceRegion[startFaceI] = -2;
1173  }
1174  else
1175  {
1176  nRegions++;
1177  }
1178  }
1179 
1180 
1181  // Check we're deciding the same on both sides. Since the regioning
1182  // is done based on nFacesPerEdge (which is synced) this should
1183  // indeed be the case.
1184 
1185  labelList nbrFaceRegion(faceRegion);
1186 
1188  (
1189  mesh_,
1190  nbrFaceRegion,
1191  false // no separation
1192  );
1193 
1194  labelList toNbrRegion(nRegions, -1);
1195 
1196  for
1197  (
1198  label faceI = mesh_.nInternalFaces();
1199  faceI < mesh_.nFaces();
1200  faceI++
1201  )
1202  {
1203  // Get the neighbouring region.
1204  label nbrRegion = nbrFaceRegion[faceI];
1205  label myRegion = faceRegion[faceI];
1206 
1207  if (myRegion <= -1 || nbrRegion <= -1)
1208  {
1209  if (nbrRegion != myRegion)
1210  {
1211  FatalErrorIn("removeFaces::setRefinement")
1212  << "Inconsistent face region across coupled patches."
1213  << endl
1214  << "This side has for faceI:" << faceI
1215  << " region:" << myRegion << endl
1216  << "The other side has region:" << nbrRegion
1217  << endl
1218  << "(region -1 means face is to be deleted)"
1219  << abort(FatalError);
1220  }
1221  }
1222  else if (toNbrRegion[myRegion] == -1)
1223  {
1224  // First visit of region. Store correspondence.
1225  toNbrRegion[myRegion] = nbrRegion;
1226  }
1227  else
1228  {
1229  // Second visit of this region.
1230  if (toNbrRegion[myRegion] != nbrRegion)
1231  {
1232  FatalErrorIn("removeFaces::setRefinement")
1233  << "Inconsistent face region across coupled patches."
1234  << endl
1235  << "This side has for faceI:" << faceI
1236  << " region:" << myRegion
1237  << " with coupled neighbouring regions:"
1238  << toNbrRegion[myRegion] << " and "
1239  << nbrRegion
1240  << abort(FatalError);
1241  }
1242  }
1243  }
1244  }
1245 
1246  //if (debug)
1247  //{
1248  // labelListList regionToFaces(invertOneToMany(nRegions, faceRegion));
1249  //
1250  // forAll(regionToFaces, regionI)
1251  // {
1252  // Pout<< " " << regionI << " faces:" << regionToFaces[regionI]
1253  // << endl;
1254  // }
1255  //}
1256 
1257 
1258  // Points to be removed
1259  // ~~~~~~~~~~~~~~~~~~~~
1260 
1261  labelHashSet pointsToRemove(4*faceLabels.size());
1262 
1263 
1264  // Per point count the number of unremoved edges. Store the ones that
1265  // are only used by 2 unremoved edges.
1266  {
1267  // Usage of points by non-removed edges.
1268  labelList nEdgesPerPoint(mesh_.nPoints());
1269 
1270  const labelListList& pointEdges = mesh_.pointEdges();
1271 
1272  forAll(pointEdges, pointI)
1273  {
1274  nEdgesPerPoint[pointI] = pointEdges[pointI].size();
1275  }
1276 
1277  forAllConstIter(labelHashSet, edgesToRemove, iter)
1278  {
1279  // Edge will get removed.
1280  const edge& e = mesh_.edges()[iter.key()];
1281 
1282  forAll(e, i)
1283  {
1284  nEdgesPerPoint[e[i]]--;
1285  }
1286  }
1287 
1288  // Check locally (before synchronizing) for strangeness
1289  forAll(nEdgesPerPoint, pointI)
1290  {
1291  if (nEdgesPerPoint[pointI] == 1)
1292  {
1293  FatalErrorIn("removeFaces::setRefinement")
1294  << "Problem : point would get 1 edge using it only."
1295  << " pointI:" << pointI
1296  << " coord:" << mesh_.points()[pointI]
1297  << abort(FatalError);
1298  }
1299  }
1300 
1301  // Synchronize point usage. This is to make sure that both sides remove
1302  // (or not remove) a point on the boundary at the same time.
1304  (
1305  mesh_,
1306  nEdgesPerPoint,
1307  maxEqOp<label>(),
1308  labelMin,
1309  false // no separation
1310  );
1311 
1312  forAll(nEdgesPerPoint, pointI)
1313  {
1314  if (nEdgesPerPoint[pointI] == 0)
1315  {
1316  pointsToRemove.insert(pointI);
1317  }
1318  else if (nEdgesPerPoint[pointI] == 1)
1319  {
1320  // Already checked before
1321  }
1322  else if (nEdgesPerPoint[pointI] == 2)
1323  {
1324  // Remove point and merge edges.
1325  pointsToRemove.insert(pointI);
1326  }
1327  }
1328  }
1329 
1330 
1331  if (debug)
1332  {
1333  OFstream str(mesh_.time().path()/"pointsToRemove.obj");
1334  Pout<< "Dumping pointsToRemove to " << str.name() << endl;
1335 
1336  forAllConstIter(labelHashSet, pointsToRemove, iter)
1337  {
1338  meshTools::writeOBJ(str, mesh_.points()[iter.key()]);
1339  }
1340  }
1341 
1342 
1343  // All faces affected in any way
1344  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1345 
1346  // Get all faces affected in any way by removal of points/edges/faces/cells
1347  boolList affectedFace
1348  (
1349  getFacesAffected
1350  (
1351  cellRegion,
1352  cellRegionMaster,
1353  faceLabels,
1354  edgesToRemove,
1355  pointsToRemove
1356  )
1357  );
1358 
1359  //
1360  // Now we know
1361  // - faceLabels : faces to remove (sync since no boundary faces)
1362  // - cellRegion/Master : cells to remove (sync since cells)
1363  // - pointsToRemove : points to remove (sync)
1364  // - faceRegion : connected face region of faces to be merged (sync)
1365  // - affectedFace : faces with points removed and/or owner/neighbour
1366  // changed (non sync)
1367 
1368 
1369  // Start modifying mesh and keep track of faces changed.
1370 
1371 
1372  // Do all removals
1373  // ~~~~~~~~~~~~~~~
1374 
1375  // Remove split faces.
1376  forAll(faceLabels, labelI)
1377  {
1378  label faceI = faceLabels[labelI];
1379 
1380  // Remove face if not yet uptodate (which is never; but want to be
1381  // consistent with rest of face removals/modifications)
1382  if (affectedFace[faceI])
1383  {
1384  affectedFace[faceI] = false;
1385 
1386  meshMod.setAction(polyRemoveFace(faceI, -1));
1387  }
1388  }
1389 
1390 
1391  // Remove points.
1392  forAllConstIter(labelHashSet, pointsToRemove, iter)
1393  {
1394  label pointI = iter.key();
1395 
1396  meshMod.setAction(polyRemovePoint(pointI, -1));
1397  }
1398 
1399 
1400  // Remove cells.
1401  forAll(cellRegion, cellI)
1402  {
1403  label region = cellRegion[cellI];
1404 
1405  if (region != -1 && (cellI != cellRegionMaster[region]))
1406  {
1407  meshMod.setAction(polyRemoveCell(cellI, cellRegionMaster[region]));
1408  }
1409  }
1410 
1411 
1412 
1413  // Merge faces across edges to be merged
1414  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1415 
1416  // Invert faceRegion so we get region to faces.
1417  {
1418  labelListList regionToFaces(invertOneToMany(nRegions, faceRegion));
1419 
1420  forAll(regionToFaces, regionI)
1421  {
1422  const labelList& rFaces = regionToFaces[regionI];
1423 
1424  if (rFaces.size() <= 1)
1425  {
1426  FatalErrorIn("setRefinement")
1427  << "Region:" << regionI
1428  << " contains only faces " << rFaces
1429  << abort(FatalError);
1430  }
1431 
1432  // rFaces[0] is master, rest gets removed.
1433  mergeFaces
1434  (
1435  cellRegion,
1436  cellRegionMaster,
1437  pointsToRemove,
1438  rFaces,
1439  meshMod
1440  );
1441 
1442  forAll(rFaces, i)
1443  {
1444  affectedFace[rFaces[i]] = false;
1445  }
1446  }
1447  }
1448 
1449 
1450  // Remaining affected faces
1451  // ~~~~~~~~~~~~~~~~~~~~~~~~
1452 
1453  // Check if any remaining faces have not been updated for new slave/master
1454  // or points removed.
1455  forAll(affectedFace, faceI)
1456  {
1457  if (affectedFace[faceI])
1458  {
1459  affectedFace[faceI] = false;
1460 
1461  face f(filterFace(pointsToRemove, faceI));
1462 
1463  label own = mesh_.faceOwner()[faceI];
1464 
1465  if (cellRegion[own] != -1)
1466  {
1467  own = cellRegionMaster[cellRegion[own]];
1468  }
1469 
1470  label patchID, zoneID, zoneFlip;
1471 
1472  getFaceInfo(faceI, patchID, zoneID, zoneFlip);
1473 
1474  label nei = -1;
1475 
1476  if (mesh_.isInternalFace(faceI))
1477  {
1478  nei = mesh_.faceNeighbour()[faceI];
1479 
1480  if (cellRegion[nei] != -1)
1481  {
1482  nei = cellRegionMaster[cellRegion[nei]];
1483  }
1484  }
1485 
1486 // if (debug)
1487 // {
1488 // Pout<< "Modifying " << faceI
1489 // << " old verts:" << mesh_.faces()[faceI]
1490 // << " for new verts:" << f
1491 // << " or for new owner " << own << " or for new nei "
1492 // << nei
1493 // << endl;
1494 // }
1495 
1496  modFace
1497  (
1498  f, // modified face
1499  faceI, // label of face being modified
1500  own, // owner
1501  nei, // neighbour
1502  false, // face flip
1503  patchID, // patch for face
1504  false, // remove from zone
1505  zoneID, // zone for face
1506  zoneFlip, // face flip in zone
1507 
1508  meshMod
1509  );
1510  }
1511  }
1512 }
1513 
1514 
1515 // ************************ vim: set sw=4 sts=4 et: ************************ //