1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "removeFaces.H"
29 #include "polyTopoChange.H"
30 #include "meshTools.H"
31 #include "polyModifyFace.H"
32 #include "polyRemoveFace.H"
33 #include "polyRemoveCell.H"
34 #include "polyRemovePoint.H"
35 #include "syncTools.H"
37 #include "indirectPrimitivePatch.H"
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
46 defineTypeNameAndDebug(removeFaces, 0);
50 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
52 // Changes region of connected set of cells. Can be recursive since hopefully
53 // only small area of faces removed in one go.
54 void Foam::removeFaces::changeCellRegion
57 const label oldRegion,
58 const label newRegion,
62 if (cellRegion[cellI] == oldRegion)
64 cellRegion[cellI] = newRegion;
66 // Step to neighbouring cells
68 const labelList& cCells = mesh_.cellCells()[cellI];
72 changeCellRegion(cCells[i], oldRegion, newRegion, cellRegion);
78 // Changes region of connected set of faces. Returns number of changed faces.
79 Foam::label Foam::removeFaces::changeFaceRegion
81 const labelList& cellRegion,
82 const boolList& removedFace,
83 const labelList& nFacesPerEdge,
85 const label newRegion,
86 const labelList& fEdges,
92 if (faceRegion[faceI] == -1 && !removedFace[faceI])
94 faceRegion[faceI] = newRegion;
98 // Storage for on-the-fly addressing
99 DynamicList<label> fe;
100 DynamicList<label> ef;
102 // Step to neighbouring faces across edges that will get removed
105 label edgeI = fEdges[i];
107 if (nFacesPerEdge[edgeI] >= 0 && nFacesPerEdge[edgeI] <= 2)
109 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
113 label nbrFaceI = eFaces[j];
115 const labelList& fEdges1 = mesh_.faceEdges(nbrFaceI, fe);
117 nChanged += changeFaceRegion
135 // Mark all faces affected in any way by
136 // - removal of cells
137 // - removal of faces
138 // - removal of edges
139 // - removal of points
140 Foam::boolList Foam::removeFaces::getFacesAffected
142 const labelList& cellRegion,
143 const labelList& cellRegionMaster,
144 const labelList& facesToRemove,
145 const labelHashSet& edgesToRemove,
146 const labelHashSet& pointsToRemove
149 boolList affectedFace(mesh_.nFaces(), false);
151 // Mark faces affected by removal of cells
152 forAll(cellRegion, cellI)
154 label region = cellRegion[cellI];
156 if (region != -1 && (cellI != cellRegionMaster[region]))
158 const labelList& cFaces = mesh_.cells()[cellI];
160 forAll(cFaces, cFaceI)
162 affectedFace[cFaces[cFaceI]] = true;
167 // Mark faces affected by removal of face.
168 forAll(facesToRemove, i)
170 affectedFace[facesToRemove[i]] = true;
173 // Mark faces affected by removal of edges
174 forAllConstIter(labelHashSet, edgesToRemove, iter)
176 const labelList& eFaces = mesh_.edgeFaces(iter.key());
178 forAll(eFaces, eFaceI)
180 affectedFace[eFaces[eFaceI]] = true;
184 // Mark faces affected by removal of points
185 forAllConstIter(labelHashSet, pointsToRemove, iter)
187 label pointI = iter.key();
189 const labelList& pFaces = mesh_.pointFaces()[pointI];
191 forAll(pFaces, pFaceI)
193 affectedFace[pFaces[pFaceI]] = true;
200 void Foam::removeFaces::writeOBJ
202 const indirectPrimitivePatch& fp,
203 const fileName& fName
207 Pout<< "removeFaces::writeOBJ : Writing faces to file "
208 << str.name() << endl;
210 const pointField& localPoints = fp.localPoints();
212 forAll(localPoints, i)
214 meshTools::writeOBJ(str, localPoints[i]);
217 const faceList& localFaces = fp.localFaces();
219 forAll(localFaces, i)
221 const face& f = localFaces[i];
227 str<< ' ' << f[fp]+1;
234 // Inserts commands to merge faceLabels into one face.
235 void Foam::removeFaces::mergeFaces
237 const labelList& cellRegion,
238 const labelList& cellRegionMaster,
239 const labelHashSet& pointsToRemove,
240 const labelList& faceLabels,
241 polyTopoChange& meshMod
244 // Construct addressing engine from faceLabels (in order of faceLabels as
246 indirectPrimitivePatch fp
256 // Get outside vertices (in local vertex numbering)
258 if (fp.edgeLoops().size() != 1)
260 writeOBJ(fp, mesh_.time().path()/"facesToBeMerged.obj");
261 FatalErrorIn("removeFaces::mergeFaces")
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);
268 const labelList& edgeLoop = fp.edgeLoops()[0];
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
275 label masterIndex = -1;
276 bool reverseLoop = false;
278 const labelList& pFaces = fp.pointFaces()[edgeLoop[0]];
280 // Find face among pFaces which uses edgeLoop[1]
283 label faceI = pFaces[i];
285 const face& f = fp.localFaces()[faceI];
287 label index1 = findIndex(f, edgeLoop[1]);
291 // Check whether consecutive to edgeLoop[0]
292 label index0 = findIndex(f, edgeLoop[0]);
296 if (index1 == f.fcIndex(index0))
302 else if (index1 == f.rcIndex(index0))
312 if (masterIndex == -1)
314 writeOBJ(fp, mesh_.time().path()/"facesToBeMerged.obj");
315 FatalErrorIn("removeFaces::mergeFaces")
316 << "Problem" << abort(FatalError);
320 // Modify the master face.
321 // ~~~~~~~~~~~~~~~~~~~~~~~
323 // Modify first face.
324 label faceI = faceLabels[masterIndex];
326 label own = mesh_.faceOwner()[faceI];
328 if (cellRegion[own] != -1)
330 own = cellRegionMaster[cellRegion[own]];
333 label patchID, zoneID, zoneFlip;
335 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
339 if (mesh_.isInternalFace(faceI))
341 nei = mesh_.faceNeighbour()[faceI];
343 if (cellRegion[nei] != -1)
345 nei = cellRegionMaster[cellRegion[nei]];
350 DynamicList<label> faceVerts(edgeLoop.size());
354 label pointI = fp.meshPoints()[edgeLoop[i]];
356 if (pointsToRemove.found(pointI))
358 //Pout<< "**Removing point " << pointI << " from "
359 // << edgeLoop << endl;
363 faceVerts.append(pointI);
368 mergedFace.transfer(faceVerts);
376 // Pout<< "Modifying masterface " << faceI
377 // << " from faces:" << faceLabels
378 // << " old verts:" << UIndirectList<face>(mesh_.faces(), faceLabels)
379 // << " for new verts:"
381 // << " possibly new owner " << own
382 // << " or new nei " << nei
388 mergedFace, // modified face
389 faceI, // label of face being modified
393 patchID, // patch for face
394 false, // remove from zone
395 zoneID, // zone for face
396 zoneFlip, // face flip in zone
402 // Remove all but master face.
403 forAll(faceLabels, patchFaceI)
405 if (patchFaceI != masterIndex)
407 //Pout<< "Removing face " << faceLabels[patchFaceI] << endl;
409 meshMod.setAction(polyRemoveFace(faceLabels[patchFaceI], faceI));
415 // Get patch, zone info for faceI
416 void Foam::removeFaces::getFaceInfo
427 if (!mesh_.isInternalFace(faceI))
429 patchID = mesh_.boundaryMesh().whichPatch(faceI);
432 zoneID = mesh_.faceZones().whichZone(faceI);
438 const faceZone& fZone = mesh_.faceZones()[zoneID];
440 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
445 // Return face with all pointsToRemove removed.
446 Foam::face Foam::removeFaces::filterFace
448 const labelHashSet& pointsToRemove,
452 const face& f = mesh_.faces()[faceI];
454 labelList newFace(f.size(), -1);
462 if (!pointsToRemove.found(vertI))
464 newFace[newFp++] = vertI;
468 newFace.setSize(newFp);
470 return face(newFace);
474 // Wrapper for meshMod.modifyFace. Reverses face if own>nei.
475 void Foam::removeFaces::modFace
478 const label masterFaceID,
481 const bool flipFaceFlux,
482 const label newPatchID,
483 const bool removeFromZone,
487 polyTopoChange& meshMod
490 if ((nei == -1) || (own < nei))
494 // Pout<< "ModifyFace (unreversed) :"
495 // << " faceI:" << masterFaceID
499 // << " flipFaceFlux:" << flipFaceFlux
500 // << " newPatchID:" << newPatchID
501 // << " removeFromZone:" << removeFromZone
502 // << " zoneID:" << zoneID
503 // << " zoneFlip:" << zoneFlip
512 masterFaceID, // label of face being modified
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
527 // Pout<< "ModifyFace (!reversed) :"
528 // << " faceI:" << masterFaceID
529 // << " f:" << f.reverseFace()
532 // << " flipFaceFlux:" << flipFaceFlux
533 // << " newPatchID:" << newPatchID
534 // << " removeFromZone:" << removeFromZone
535 // << " zoneID:" << zoneID
536 // << " zoneFlip:" << zoneFlip
544 f.reverseFace(),// modified face
545 masterFaceID, // label of face being modified
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
559 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
561 // Construct from mesh
562 Foam::removeFaces::removeFaces
564 const polyMesh& mesh,
573 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
575 // Removing face connects cells. This function works out a consistent set of
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.
581 Foam::label Foam::removeFaces::compatibleRemoves
583 const labelList& facesToRemove,
584 labelList& cellRegion,
585 labelList& regionMaster,
586 labelList& newFacesToRemove
589 const labelList& faceOwner = mesh_.faceOwner();
590 const labelList& faceNeighbour = mesh_.faceNeighbour();
592 cellRegion.setSize(mesh_.nCells());
595 regionMaster.setSize(mesh_.nCells());
600 forAll(facesToRemove, i)
602 label faceI = facesToRemove[i];
604 if (!mesh_.isInternalFace(faceI))
608 "removeFaces::compatibleRemoves(const labelList&"
609 ", labelList&, labelList&, labelList&)"
610 ) << "Not internal face:" << faceI << abort(FatalError);
614 label own = faceOwner[faceI];
615 label nei = faceNeighbour[faceI];
617 label region0 = cellRegion[own];
618 label region1 = cellRegion[nei];
625 cellRegion[own] = nRegions;
626 cellRegion[nei] = nRegions;
628 // Make owner (lowest numbered!) the master of the region
629 regionMaster[nRegions] = own;
634 // Add owner to neighbour region
635 cellRegion[own] = region1;
636 // See if owner becomes the master of the region
637 regionMaster[region1] = min(own, regionMaster[region1]);
644 // Add neighbour to owner region
645 cellRegion[nei] = region0;
646 // nei is higher numbered than own so guaranteed not lower
647 // than master of region0.
649 else if (region0 != region1)
651 // Both have regions. Keep lowest numbered region and master.
652 label freedRegion = -1;
653 label keptRegion = -1;
655 if (region0 < region1)
660 region1, // old region
661 region0, // new region
665 keptRegion = region0;
666 freedRegion = region1;
668 else if (region1 < region0)
673 region0, // old region
674 region1, // new region
678 keptRegion = region1;
679 freedRegion = region0;
682 label master0 = regionMaster[region0];
683 label master1 = regionMaster[region1];
685 regionMaster[freedRegion] = -1;
686 regionMaster[keptRegion] = min(master0, master1);
691 regionMaster.setSize(nRegions);
695 // - master is lowest numbered in any region
696 // - regions have more than 1 cell
698 labelList nCells(regionMaster.size(), 0);
700 forAll(cellRegion, cellI)
702 label r = cellRegion[cellI];
708 if (cellI < regionMaster[r])
712 "removeFaces::compatibleRemoves(const labelList&"
713 ", labelList&, labelList&, labelList&)"
714 ) << "Not lowest numbered : cell:" << cellI
716 << " regionmaster:" << regionMaster[r]
717 << abort(FatalError);
722 forAll(nCells, region)
724 if (nCells[region] == 1)
728 "removeFaces::compatibleRemoves(const labelList&"
729 ", labelList&, labelList&, labelList&)"
730 ) << "Region " << region
731 << " has only " << nCells[region] << " cells in it"
732 << abort(FatalError);
738 // Count number of used regions
739 label nUsedRegions = 0;
741 forAll(regionMaster, i)
743 if (regionMaster[i] != -1)
749 // Recreate facesToRemove to be consistent with the cellRegions.
750 DynamicList<label> allFacesToRemove(facesToRemove.size());
752 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
754 label own = faceOwner[faceI];
755 label nei = faceNeighbour[faceI];
757 if (cellRegion[own] != -1 && cellRegion[own] == cellRegion[nei])
759 // Both will become the same cell so add face to list of faces
761 allFacesToRemove.append(faceI);
765 newFacesToRemove.transfer(allFacesToRemove);
771 void Foam::removeFaces::setRefinement
773 const labelList& faceLabels,
774 const labelList& cellRegion,
775 const labelList& cellRegionMaster,
776 polyTopoChange& meshMod
781 faceSet facesToRemove(mesh_, "facesToRemove", faceLabels);
782 Pout<< "Writing faces to remove to faceSet " << facesToRemove.name()
784 facesToRemove.write();
787 // Make map of all faces to be removed
788 boolList removedFace(mesh_.nFaces(), false);
790 forAll(faceLabels, i)
792 label faceI = faceLabels[i];
794 if (!mesh_.isInternalFace(faceI))
798 "removeFaces::setRefinement(const labelList&"
799 ", const labelList&, const labelList&, polyTopoChange&)"
800 ) << "Face to remove is not internal face:" << faceI
801 << abort(FatalError);
804 removedFace[faceI] = true;
808 // Edges to be removed
809 // ~~~~~~~~~~~~~~~~~~~
813 labelHashSet edgesToRemove(faceLabels.size());
815 // Per face the region it is in. -1 for removed faces, -2 for regions
816 // consisting of single face only.
817 labelList faceRegion(mesh_.nFaces(), -1);
819 // Number of connected face regions
822 // Storage for on-the-fly addressing
823 DynamicList<label> fe;
824 DynamicList<label> ef;
828 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
830 // Usage of edges by non-removed faces.
831 // See below about initialization.
832 labelList nFacesPerEdge(mesh_.nEdges(), -1);
834 // Count usage of edges by non-removed faces.
835 forAll(faceLabels, i)
837 label faceI = faceLabels[i];
839 const labelList& fEdges = mesh_.faceEdges(faceI, fe);
843 label edgeI = fEdges[i];
845 if (nFacesPerEdge[edgeI] == -1)
847 nFacesPerEdge[edgeI] = mesh_.edgeFaces(edgeI, ef).size()-1;
851 nFacesPerEdge[edgeI]--;
856 // Count usage for edges not on faces-to-be-removed.
857 // Note that this only needs to be done for possibly coupled edges
858 // so we could choose to loop only over boundary faces and use faceEdges
861 forAll(mesh_.edges(), edgeI)
863 if (nFacesPerEdge[edgeI] == -1)
865 // Edge not yet handled in loop above so is not used by any
866 // face to be removed.
868 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
870 if (eFaces.size() > 2)
872 nFacesPerEdge[edgeI] = eFaces.size();
874 else if (eFaces.size() == 2)
876 // nFacesPerEdge already -1 so do nothing.
880 const edge& e = mesh_.edges()[edgeI];
882 FatalErrorIn("removeFaces::setRefinement")
883 << "Problem : edge has too few face neighbours:"
887 << " coords:" << mesh_.points()[e[0]]
888 << mesh_.points()[e[1]]
889 << abort(FatalError);
898 OFstream str(mesh_.time().path()/"edgesWithTwoFaces.obj");
899 Pout<< "Dumping edgesWithTwoFaces to " << str.name() << endl;
902 forAll(nFacesPerEdge, edgeI)
904 if (nFacesPerEdge[edgeI] == 2)
906 // Edge will get removed.
907 const edge& e = mesh_.edges()[edgeI];
909 meshTools::writeOBJ(str, mesh_.points()[e[0]]);
911 meshTools::writeOBJ(str, mesh_.points()[e[1]]);
913 str<< "l " << vertI-1 << ' ' << vertI << nl;
919 // Now all unaffected edges will have labelMax, all affected edges the
920 // number of unremoved faces.
922 // Filter for edges inbetween two remaining boundary faces that
923 // make too big an angle.
924 forAll(nFacesPerEdge, edgeI)
926 if (nFacesPerEdge[edgeI] == 2)
928 // See if they are two boundary faces
932 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
936 label faceI = eFaces[i];
938 if (!removedFace[faceI] && !mesh_.isInternalFace(faceI))
952 if (f0 != -1 && f1 != -1)
954 // Edge has two boundary faces remaining.
955 // See if should be merged.
957 label patch0 = patches.whichPatch(f0);
958 label patch1 = patches.whichPatch(f1);
960 if (patch0 != patch1)
962 // Different patches. Do not merge edge.
963 WarningIn("removeFaces::setRefinement")
964 << "not merging faces " << f0 << " and "
965 << f1 << " across patch boundary edge " << edgeI
968 // Mark so it gets preserved
969 nFacesPerEdge[edgeI] = 3;
971 else if (minCos_ < 1 && minCos_ > -1)
973 const polyPatch& pp0 = patches[patch0];
974 const vectorField& n0 = pp0.faceNormals();
981 & n0[f1 - pp0.start()]
986 WarningIn("removeFaces::setRefinement")
987 << "not merging faces " << f0 << " and "
988 << f1 << " across edge " << edgeI
991 // Angle between two remaining faces too large.
992 // Mark so it gets preserved
993 nFacesPerEdge[edgeI] = 3;
997 else if (f0 != -1 || f1 != -1)
999 const edge& e = mesh_.edges()[edgeI];
1001 // Only found one boundary face. Problem.
1002 FatalErrorIn("removeFaces::setRefinement")
1003 << "Problem : edge would have one boundary face"
1004 << " and one internal face using it." << endl
1005 << "Your remove pattern is probably incorrect." << endl
1007 << " nFaces:" << nFacesPerEdge[edgeI]
1008 << " vertices:" << e
1009 << " coords:" << mesh_.points()[e[0]]
1010 << mesh_.points()[e[1]]
1013 << abort(FatalError);
1020 // Check locally (before synchronizing) for strangeness
1021 forAll(nFacesPerEdge, edgeI)
1023 if (nFacesPerEdge[edgeI] == 1)
1025 const edge& e = mesh_.edges()[edgeI];
1027 FatalErrorIn("removeFaces::setRefinement")
1028 << "Problem : edge would get 1 face using it only"
1029 << " edge:" << edgeI
1030 << " nFaces:" << nFacesPerEdge[edgeI]
1031 << " vertices:" << e
1032 << " coords:" << mesh_.points()[e[0]]
1033 << ' ' << mesh_.points()[e[1]]
1034 << abort(FatalError);
1036 // Could check here for boundary edge with <=1 faces remaining.
1040 // Synchronize edge usage. This is to make sure that both sides remove
1041 // (or not remove) an edge on the boundary at the same time.
1043 // Coupled edges (edge0, edge1 are opposite each other)
1044 // a. edge not on face to be removed, edge has >= 3 faces
1045 // b. ,, edge has 2 faces
1046 // c. edge has >= 3 remaining faces
1047 // d. edge has 2 remaining faces (assume angle>minCos already handled)
1049 // - a + a: do not remove edge
1050 // - a + b: do not remove edge
1051 // - a + c: do not remove edge
1052 // - a + d: do not remove edge
1054 // - b + b: do not remove edge
1055 // - b + c: do not remove edge
1056 // - b + d: remove edge
1058 // - c + c: do not remove edge
1059 // - c + d: do not remove edge
1060 // - d + d: remove edge
1063 // So code situation a. with >= 3
1067 // then do max and check result.
1069 // a+a : max(3,3) = 3. do not remove
1070 // a+b : max(3,-1) = 3. do not remove
1071 // a+c : max(3,3) = 3. do not remove
1072 // a+d : max(3,2) = 3. do not remove
1074 // b+b : max(-1,-1) = -1. do not remove
1075 // b+c : max(-1,3) = 3. do not remove
1076 // b+d : max(-1,2) = 2. remove
1078 // c+c : max(3,3) = 3. do not remove
1079 // c+d : max(3,2) = 3. do not remove
1081 // d+d : max(2,2) = 2. remove
1083 syncTools::syncEdgeList
1088 labelMin, // guaranteed to be overridden by maxEqOp
1089 false // no separation
1092 // Convert to labelHashSet
1093 forAll(nFacesPerEdge, edgeI)
1095 if (nFacesPerEdge[edgeI] == 0)
1097 // 0: edge not used anymore.
1098 edgesToRemove.insert(edgeI);
1100 else if (nFacesPerEdge[edgeI] == 1)
1102 // 1: illegal. Tested above.
1104 else if (nFacesPerEdge[edgeI] == 2)
1107 edgesToRemove.insert(edgeI);
1113 OFstream str(mesh_.time().path()/"edgesToRemove.obj");
1114 Pout<< "Dumping edgesToRemove to " << str.name() << endl;
1117 forAllConstIter(labelHashSet, edgesToRemove, iter)
1119 // Edge will get removed.
1120 const edge& e = mesh_.edges()[iter.key()];
1122 meshTools::writeOBJ(str, mesh_.points()[e[0]]);
1124 meshTools::writeOBJ(str, mesh_.points()[e[1]]);
1126 str<< "l " << vertI-1 << ' ' << vertI << nl;
1131 // Walk to fill faceRegion with faces that will be connected across
1132 // edges that will be removed.
1134 label startFaceI = 0;
1138 // Find unset region.
1139 for (; startFaceI < mesh_.nFaces(); startFaceI++)
1141 if (faceRegion[startFaceI] == -1 && !removedFace[startFaceI])
1147 if (startFaceI == mesh_.nFaces())
1152 // Start walking face-edge-face, crossing edges that will get
1153 // removed. Every thus connected region will get single region
1155 label nRegion = changeFaceRegion
1162 mesh_.faceEdges(startFaceI, fe),
1168 FatalErrorIn("setRefinement") << "Problem" << abort(FatalError);
1170 else if (nRegion == 1)
1172 // Reset face to be single region.
1173 faceRegion[startFaceI] = -2;
1182 // Check we're deciding the same on both sides. Since the regioning
1183 // is done based on nFacesPerEdge (which is synced) this should
1184 // indeed be the case.
1186 labelList nbrFaceRegion(faceRegion);
1188 syncTools::swapFaceList
1192 false // no separation
1195 labelList toNbrRegion(nRegions, -1);
1199 label faceI = mesh_.nInternalFaces();
1200 faceI < mesh_.nFaces();
1204 // Get the neighbouring region.
1205 label nbrRegion = nbrFaceRegion[faceI];
1206 label myRegion = faceRegion[faceI];
1208 if (myRegion <= -1 || nbrRegion <= -1)
1210 if (nbrRegion != myRegion)
1212 FatalErrorIn("removeFaces::setRefinement")
1213 << "Inconsistent face region across coupled patches."
1215 << "This side has for faceI:" << faceI
1216 << " region:" << myRegion << endl
1217 << "The other side has region:" << nbrRegion
1219 << "(region -1 means face is to be deleted)"
1220 << abort(FatalError);
1223 else if (toNbrRegion[myRegion] == -1)
1225 // First visit of region. Store correspondence.
1226 toNbrRegion[myRegion] = nbrRegion;
1230 // Second visit of this region.
1231 if (toNbrRegion[myRegion] != nbrRegion)
1233 FatalErrorIn("removeFaces::setRefinement")
1234 << "Inconsistent face region across coupled patches."
1236 << "This side has for faceI:" << faceI
1237 << " region:" << myRegion
1238 << " with coupled neighbouring regions:"
1239 << toNbrRegion[myRegion] << " and "
1241 << abort(FatalError);
1249 // labelListList regionToFaces(invertOneToMany(nRegions, faceRegion));
1251 // forAll(regionToFaces, regionI)
1253 // Pout<< " " << regionI << " faces:" << regionToFaces[regionI]
1259 // Points to be removed
1260 // ~~~~~~~~~~~~~~~~~~~~
1262 labelHashSet pointsToRemove(4*faceLabels.size());
1265 // Per point count the number of unremoved edges. Store the ones that
1266 // are only used by 2 unremoved edges.
1268 // Usage of points by non-removed edges.
1269 labelList nEdgesPerPoint(mesh_.nPoints());
1271 const labelListList& pointEdges = mesh_.pointEdges();
1273 forAll(pointEdges, pointI)
1275 nEdgesPerPoint[pointI] = pointEdges[pointI].size();
1278 forAllConstIter(labelHashSet, edgesToRemove, iter)
1280 // Edge will get removed.
1281 const edge& e = mesh_.edges()[iter.key()];
1285 nEdgesPerPoint[e[i]]--;
1289 // Check locally (before synchronizing) for strangeness
1290 forAll(nEdgesPerPoint, pointI)
1292 if (nEdgesPerPoint[pointI] == 1)
1294 FatalErrorIn("removeFaces::setRefinement")
1295 << "Problem : point would get 1 edge using it only."
1296 << " pointI:" << pointI
1297 << " coord:" << mesh_.points()[pointI]
1298 << abort(FatalError);
1302 // Synchronize point usage. This is to make sure that both sides remove
1303 // (or not remove) a point on the boundary at the same time.
1304 syncTools::syncPointList
1310 false // no separation
1313 forAll(nEdgesPerPoint, pointI)
1315 if (nEdgesPerPoint[pointI] == 0)
1317 pointsToRemove.insert(pointI);
1319 else if (nEdgesPerPoint[pointI] == 1)
1321 // Already checked before
1323 else if (nEdgesPerPoint[pointI] == 2)
1325 // Remove point and merge edges.
1326 pointsToRemove.insert(pointI);
1334 OFstream str(mesh_.time().path()/"pointsToRemove.obj");
1335 Pout<< "Dumping pointsToRemove to " << str.name() << endl;
1337 forAllConstIter(labelHashSet, pointsToRemove, iter)
1339 meshTools::writeOBJ(str, mesh_.points()[iter.key()]);
1344 // All faces affected in any way
1345 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1347 // Get all faces affected in any way by removal of points/edges/faces/cells
1348 boolList affectedFace
1362 // - faceLabels : faces to remove (sync since no boundary faces)
1363 // - cellRegion/Master : cells to remove (sync since cells)
1364 // - pointsToRemove : points to remove (sync)
1365 // - faceRegion : connected face region of faces to be merged (sync)
1366 // - affectedFace : faces with points removed and/or owner/neighbour
1367 // changed (non sync)
1370 // Start modifying mesh and keep track of faces changed.
1376 // Remove split faces.
1377 forAll(faceLabels, labelI)
1379 label faceI = faceLabels[labelI];
1381 // Remove face if not yet uptodate (which is never; but want to be
1382 // consistent with rest of face removals/modifications)
1383 if (affectedFace[faceI])
1385 affectedFace[faceI] = false;
1387 meshMod.setAction(polyRemoveFace(faceI, -1));
1393 forAllConstIter(labelHashSet, pointsToRemove, iter)
1395 label pointI = iter.key();
1397 meshMod.setAction(polyRemovePoint(pointI, -1));
1402 forAll(cellRegion, cellI)
1404 label region = cellRegion[cellI];
1406 if (region != -1 && (cellI != cellRegionMaster[region]))
1408 meshMod.setAction(polyRemoveCell(cellI, cellRegionMaster[region]));
1414 // Merge faces across edges to be merged
1415 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1417 // Invert faceRegion so we get region to faces.
1419 labelListList regionToFaces(invertOneToMany(nRegions, faceRegion));
1421 forAll(regionToFaces, regionI)
1423 const labelList& rFaces = regionToFaces[regionI];
1425 if (rFaces.size() <= 1)
1427 FatalErrorIn("setRefinement")
1428 << "Region:" << regionI
1429 << " contains only faces " << rFaces
1430 << abort(FatalError);
1433 // rFaces[0] is master, rest gets removed.
1445 affectedFace[rFaces[i]] = false;
1451 // Remaining affected faces
1452 // ~~~~~~~~~~~~~~~~~~~~~~~~
1454 // Check if any remaining faces have not been updated for new slave/master
1455 // or points removed.
1456 forAll(affectedFace, faceI)
1458 if (affectedFace[faceI])
1460 affectedFace[faceI] = false;
1462 face f(filterFace(pointsToRemove, faceI));
1464 label own = mesh_.faceOwner()[faceI];
1466 if (cellRegion[own] != -1)
1468 own = cellRegionMaster[cellRegion[own]];
1471 label patchID, zoneID, zoneFlip;
1473 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
1477 if (mesh_.isInternalFace(faceI))
1479 nei = mesh_.faceNeighbour()[faceI];
1481 if (cellRegion[nei] != -1)
1483 nei = cellRegionMaster[cellRegion[nei]];
1489 // Pout<< "Modifying " << faceI
1490 // << " old verts:" << mesh_.faces()[faceI]
1491 // << " for new verts:" << f
1492 // << " or for new owner " << own << " or for new nei "
1500 faceI, // label of face being modified
1504 patchID, // patch for face
1505 false, // remove from zone
1506 zoneID, // zone for face
1507 zoneFlip, // face flip in zone
1516 // ************************************************************************* //