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) 2011 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 "removeFaces.H"
27 #include "polyMesh.H"
28 #include "polyTopoChange.H"
29 #include "meshTools.H"
30 #include "polyModifyFace.H"
31 #include "polyRemoveFace.H"
32 #include "polyRemoveCell.H"
33 #include "polyRemovePoint.H"
34 #include "syncTools.H"
35 #include "OFstream.H"
36 #include "indirectPrimitivePatch.H"
37 #include "Time.H"
38 #include "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 
199 void Foam::removeFaces::writeOBJ
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  );
1089 
1090  // Convert to labelHashSet
1091  forAll(nFacesPerEdge, edgeI)
1092  {
1093  if (nFacesPerEdge[edgeI] == 0)
1094  {
1095  // 0: edge not used anymore.
1096  edgesToRemove.insert(edgeI);
1097  }
1098  else if (nFacesPerEdge[edgeI] == 1)
1099  {
1100  // 1: illegal. Tested above.
1101  }
1102  else if (nFacesPerEdge[edgeI] == 2)
1103  {
1104  // 2: merge faces.
1105  edgesToRemove.insert(edgeI);
1106  }
1107  }
1108 
1109  if (debug)
1110  {
1111  OFstream str(mesh_.time().path()/"edgesToRemove.obj");
1112  Pout<< "Dumping edgesToRemove to " << str.name() << endl;
1113  label vertI = 0;
1114 
1115  forAllConstIter(labelHashSet, edgesToRemove, iter)
1116  {
1117  // Edge will get removed.
1118  const edge& e = mesh_.edges()[iter.key()];
1119 
1120  meshTools::writeOBJ(str, mesh_.points()[e[0]]);
1121  vertI++;
1122  meshTools::writeOBJ(str, mesh_.points()[e[1]]);
1123  vertI++;
1124  str<< "l " << vertI-1 << ' ' << vertI << nl;
1125  }
1126  }
1127 
1128 
1129  // Walk to fill faceRegion with faces that will be connected across
1130  // edges that will be removed.
1131 
1132  label startFaceI = 0;
1133 
1134  while (true)
1135  {
1136  // Find unset region.
1137  for (; startFaceI < mesh_.nFaces(); startFaceI++)
1138  {
1139  if (faceRegion[startFaceI] == -1 && !removedFace[startFaceI])
1140  {
1141  break;
1142  }
1143  }
1144 
1145  if (startFaceI == mesh_.nFaces())
1146  {
1147  break;
1148  }
1149 
1150  // Start walking face-edge-face, crossing edges that will get
1151  // removed. Every thus connected region will get single region
1152  // number.
1153  label nRegion = changeFaceRegion
1154  (
1155  cellRegion,
1156  removedFace,
1157  nFacesPerEdge,
1158  startFaceI,
1159  nRegions,
1160  mesh_.faceEdges(startFaceI, fe),
1161  faceRegion
1162  );
1163 
1164  if (nRegion < 1)
1165  {
1166  FatalErrorIn("setRefinement") << "Problem" << abort(FatalError);
1167  }
1168  else if (nRegion == 1)
1169  {
1170  // Reset face to be single region.
1171  faceRegion[startFaceI] = -2;
1172  }
1173  else
1174  {
1175  nRegions++;
1176  }
1177  }
1178 
1179 
1180  // Check we're deciding the same on both sides. Since the regioning
1181  // is done based on nFacesPerEdge (which is synced) this should
1182  // indeed be the case.
1183 
1184  labelList nbrFaceRegion(faceRegion);
1185 
1187  (
1188  mesh_,
1189  nbrFaceRegion
1190  );
1191 
1192  labelList toNbrRegion(nRegions, -1);
1193 
1194  for
1195  (
1196  label faceI = mesh_.nInternalFaces();
1197  faceI < mesh_.nFaces();
1198  faceI++
1199  )
1200  {
1201  // Get the neighbouring region.
1202  label nbrRegion = nbrFaceRegion[faceI];
1203  label myRegion = faceRegion[faceI];
1204 
1205  if (myRegion <= -1 || nbrRegion <= -1)
1206  {
1207  if (nbrRegion != myRegion)
1208  {
1209  FatalErrorIn("removeFaces::setRefinement")
1210  << "Inconsistent face region across coupled patches."
1211  << endl
1212  << "This side has for faceI:" << faceI
1213  << " region:" << myRegion << endl
1214  << "The other side has region:" << nbrRegion
1215  << endl
1216  << "(region -1 means face is to be deleted)"
1217  << abort(FatalError);
1218  }
1219  }
1220  else if (toNbrRegion[myRegion] == -1)
1221  {
1222  // First visit of region. Store correspondence.
1223  toNbrRegion[myRegion] = nbrRegion;
1224  }
1225  else
1226  {
1227  // Second visit of this region.
1228  if (toNbrRegion[myRegion] != nbrRegion)
1229  {
1230  FatalErrorIn("removeFaces::setRefinement")
1231  << "Inconsistent face region across coupled patches."
1232  << endl
1233  << "This side has for faceI:" << faceI
1234  << " region:" << myRegion
1235  << " with coupled neighbouring regions:"
1236  << toNbrRegion[myRegion] << " and "
1237  << nbrRegion
1238  << abort(FatalError);
1239  }
1240  }
1241  }
1242  }
1243 
1244  //if (debug)
1245  //{
1246  // labelListList regionToFaces(invertOneToMany(nRegions, faceRegion));
1247  //
1248  // forAll(regionToFaces, regionI)
1249  // {
1250  // Pout<< " " << regionI << " faces:" << regionToFaces[regionI]
1251  // << endl;
1252  // }
1253  //}
1254 
1255 
1256  // Points to be removed
1257  // ~~~~~~~~~~~~~~~~~~~~
1258 
1259  labelHashSet pointsToRemove(4*faceLabels.size());
1260 
1261 
1262  // Per point count the number of unremoved edges. Store the ones that
1263  // are only used by 2 unremoved edges.
1264  {
1265  // Usage of points by non-removed edges.
1266  labelList nEdgesPerPoint(mesh_.nPoints());
1267 
1268  const labelListList& pointEdges = mesh_.pointEdges();
1269 
1270  forAll(pointEdges, pointI)
1271  {
1272  nEdgesPerPoint[pointI] = pointEdges[pointI].size();
1273  }
1274 
1275  forAllConstIter(labelHashSet, edgesToRemove, iter)
1276  {
1277  // Edge will get removed.
1278  const edge& e = mesh_.edges()[iter.key()];
1279 
1280  forAll(e, i)
1281  {
1282  nEdgesPerPoint[e[i]]--;
1283  }
1284  }
1285 
1286  // Check locally (before synchronizing) for strangeness
1287  forAll(nEdgesPerPoint, pointI)
1288  {
1289  if (nEdgesPerPoint[pointI] == 1)
1290  {
1291  FatalErrorIn("removeFaces::setRefinement")
1292  << "Problem : point would get 1 edge using it only."
1293  << " pointI:" << pointI
1294  << " coord:" << mesh_.points()[pointI]
1295  << abort(FatalError);
1296  }
1297  }
1298 
1299  // Synchronize point usage. This is to make sure that both sides remove
1300  // (or not remove) a point on the boundary at the same time.
1302  (
1303  mesh_,
1304  nEdgesPerPoint,
1305  maxEqOp<label>(),
1306  labelMin
1307  );
1308 
1309  forAll(nEdgesPerPoint, pointI)
1310  {
1311  if (nEdgesPerPoint[pointI] == 0)
1312  {
1313  pointsToRemove.insert(pointI);
1314  }
1315  else if (nEdgesPerPoint[pointI] == 1)
1316  {
1317  // Already checked before
1318  }
1319  else if (nEdgesPerPoint[pointI] == 2)
1320  {
1321  // Remove point and merge edges.
1322  pointsToRemove.insert(pointI);
1323  }
1324  }
1325  }
1326 
1327 
1328  if (debug)
1329  {
1330  OFstream str(mesh_.time().path()/"pointsToRemove.obj");
1331  Pout<< "Dumping pointsToRemove to " << str.name() << endl;
1332 
1333  forAllConstIter(labelHashSet, pointsToRemove, iter)
1334  {
1335  meshTools::writeOBJ(str, mesh_.points()[iter.key()]);
1336  }
1337  }
1338 
1339 
1340  // All faces affected in any way
1341  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1342 
1343  // Get all faces affected in any way by removal of points/edges/faces/cells
1344  boolList affectedFace
1345  (
1346  getFacesAffected
1347  (
1348  cellRegion,
1349  cellRegionMaster,
1350  faceLabels,
1351  edgesToRemove,
1352  pointsToRemove
1353  )
1354  );
1355 
1356  //
1357  // Now we know
1358  // - faceLabels : faces to remove (sync since no boundary faces)
1359  // - cellRegion/Master : cells to remove (sync since cells)
1360  // - pointsToRemove : points to remove (sync)
1361  // - faceRegion : connected face region of faces to be merged (sync)
1362  // - affectedFace : faces with points removed and/or owner/neighbour
1363  // changed (non sync)
1364 
1365 
1366  // Start modifying mesh and keep track of faces changed.
1367 
1368 
1369  // Do all removals
1370  // ~~~~~~~~~~~~~~~
1371 
1372  // Remove split faces.
1373  forAll(faceLabels, labelI)
1374  {
1375  label faceI = faceLabels[labelI];
1376 
1377  // Remove face if not yet uptodate (which is never; but want to be
1378  // consistent with rest of face removals/modifications)
1379  if (affectedFace[faceI])
1380  {
1381  affectedFace[faceI] = false;
1382 
1383  meshMod.setAction(polyRemoveFace(faceI, -1));
1384  }
1385  }
1386 
1387 
1388  // Remove points.
1389  forAllConstIter(labelHashSet, pointsToRemove, iter)
1390  {
1391  label pointI = iter.key();
1392 
1393  meshMod.setAction(polyRemovePoint(pointI, -1));
1394  }
1395 
1396 
1397  // Remove cells.
1398  forAll(cellRegion, cellI)
1399  {
1400  label region = cellRegion[cellI];
1401 
1402  if (region != -1 && (cellI != cellRegionMaster[region]))
1403  {
1404  meshMod.setAction(polyRemoveCell(cellI, cellRegionMaster[region]));
1405  }
1406  }
1407 
1408 
1409 
1410  // Merge faces across edges to be merged
1411  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1412 
1413  // Invert faceRegion so we get region to faces.
1414  {
1415  labelListList regionToFaces(invertOneToMany(nRegions, faceRegion));
1416 
1417  forAll(regionToFaces, regionI)
1418  {
1419  const labelList& rFaces = regionToFaces[regionI];
1420 
1421  if (rFaces.size() <= 1)
1422  {
1423  FatalErrorIn("setRefinement")
1424  << "Region:" << regionI
1425  << " contains only faces " << rFaces
1426  << abort(FatalError);
1427  }
1428 
1429  // rFaces[0] is master, rest gets removed.
1430  mergeFaces
1431  (
1432  cellRegion,
1433  cellRegionMaster,
1434  pointsToRemove,
1435  rFaces,
1436  meshMod
1437  );
1438 
1439  forAll(rFaces, i)
1440  {
1441  affectedFace[rFaces[i]] = false;
1442  }
1443  }
1444  }
1445 
1446 
1447  // Remaining affected faces
1448  // ~~~~~~~~~~~~~~~~~~~~~~~~
1449 
1450  // Check if any remaining faces have not been updated for new slave/master
1451  // or points removed.
1452  forAll(affectedFace, faceI)
1453  {
1454  if (affectedFace[faceI])
1455  {
1456  affectedFace[faceI] = false;
1457 
1458  face f(filterFace(pointsToRemove, faceI));
1459 
1460  label own = mesh_.faceOwner()[faceI];
1461 
1462  if (cellRegion[own] != -1)
1463  {
1464  own = cellRegionMaster[cellRegion[own]];
1465  }
1466 
1467  label patchID, zoneID, zoneFlip;
1468 
1469  getFaceInfo(faceI, patchID, zoneID, zoneFlip);
1470 
1471  label nei = -1;
1472 
1473  if (mesh_.isInternalFace(faceI))
1474  {
1475  nei = mesh_.faceNeighbour()[faceI];
1476 
1477  if (cellRegion[nei] != -1)
1478  {
1479  nei = cellRegionMaster[cellRegion[nei]];
1480  }
1481  }
1482 
1483 // if (debug)
1484 // {
1485 // Pout<< "Modifying " << faceI
1486 // << " old verts:" << mesh_.faces()[faceI]
1487 // << " for new verts:" << f
1488 // << " or for new owner " << own << " or for new nei "
1489 // << nei
1490 // << endl;
1491 // }
1492 
1493  modFace
1494  (
1495  f, // modified face
1496  faceI, // label of face being modified
1497  own, // owner
1498  nei, // neighbour
1499  false, // face flip
1500  patchID, // patch for face
1501  false, // remove from zone
1502  zoneID, // zone for face
1503  zoneFlip, // face flip in zone
1504 
1505  meshMod
1506  );
1507  }
1508  }
1509 }
1510 
1511 
1512 // ************************************************************************* //
Output to file stream.
Definition: OFstream.H:81
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
A list of face labels.
Definition: faceSet.H:48
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:390
dimensioned< scalar > mag(const dimensioned< Type > &)
labelList f(nPoints)
Class containing data for point removal.
labelListList invertOneToMany(const label len, const labelUList &)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:67
Class containing data for cell removal.
scalar f1
Definition: createFields.H:28
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
List< face > faceList
Definition: faceListFwd.H:43
static const labelSphericalTensor labelI(1)
Identity labelTensor.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Namespace for OpenFOAM.
Class containing data for face removal.
label compatibleRemoves(const labelList &inPiercedFaces, labelList &cellRegion, labelList &cellRegionMaster, labelList &outPiercedFaces) const
Find faces including those with cells which have the same mastercell.
Definition: removeFaces.C:581
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
patches[0]
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
virtual bool write() const
Write using setting from DB.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
static const char nl
Definition: Ostream.H:260
const double e
Elementary charge.
Definition: doubleFloat.H:78
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define WarningIn(functionName)
Report a warning using Foam::Warning.
Foam::polyBoundaryMesh.
#define forAll(list, i)
Definition: UList.H:421
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
static const label labelMin
Definition: label.H:61
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const word & name() const
Return name.
Definition: IOobject.H:260
const Field< PointType > & faceNormals() const
Return face normals for patch.
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:88
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
void setRefinement(const labelList &piercedFaces, const labelList &cellRegion, const labelList &cellRegionMaster, polyTopoChange &) const
Play commands into polyTopoChange to remove faces.
Definition: removeFaces.C:771
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
void reverse(UList< T > &, const label n)
Definition: UListI.H:322
Direct mesh changes based on v1.3 polyTopoChange syntax.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:97
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:209
static void swapFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled face values.
Definition: syncTools.H:462
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
List< bool > boolList
Bool container classes.
Definition: boolList.H:50