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