52Foam::label Foam::meshCutAndRemove::firstCommon
54 const labelList& elems1,
55 const labelList& elems2
60 label index1 = elems2.find(elems1[elemI]);
72bool Foam::meshCutAndRemove::isIn
78 label index = cuts.find(twoCuts[0]);
87 cuts[cuts.fcIndex(index)] == twoCuts[1]
88 || cuts[cuts.rcIndex(index)] == twoCuts[1]
95Foam::label Foam::meshCutAndRemove::findCutCell
101 forAll(cellLabels, labelI)
103 label celli = cellLabels[labelI];
105 if (cuts.cellLoops()[celli].size())
114Foam::label Foam::meshCutAndRemove::findInternalFacePoint
127 label facei =
pFaces[pFacei];
129 if (
mesh().isInternalFace(facei))
146Foam::label Foam::meshCutAndRemove::findPatchFacePoint
149 const label exposedPatchi
157 label pointi =
f[fp];
178 const cellCuts& cuts,
179 const label exposedPatchi,
193 if (cellLoops[own].size() && firstCommon(
f, anchorPts[own]) == -1)
201 if (
mesh().isInternalFace(facei))
205 if (cellLoops[nei].size() && firstCommon(
f, anchorPts[nei]) == -1)
213 if (
patchID == -1 && (own == -1 || nei == -1))
221void Foam::meshCutAndRemove::getZoneInfo
236 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
241void Foam::meshCutAndRemove::addFace
243 polyTopoChange& meshMod,
245 const label masterPointi,
255 getZoneInfo(facei,
zoneID, zoneFlip);
257 if ((nei == -1) || (own != -1 && own < nei))
262 Pout<<
"Adding face " << newFace
263 <<
" with new owner:" << own
264 <<
" with new neighbour:" << nei
266 <<
" anchor:" << masterPointi
268 <<
" zoneFlip:" << zoneFlip
294 Pout<<
"Adding (reversed) face " << newFace.reverseFace()
295 <<
" with new owner:" << nei
296 <<
" with new neighbour:" << own
298 <<
" anchor:" << masterPointi
300 <<
" zoneFlip:" << zoneFlip
308 newFace.reverseFace(),
325void Foam::meshCutAndRemove::modFace
327 polyTopoChange& meshMod,
338 getZoneInfo(facei,
zoneID, zoneFlip);
342 (own !=
mesh().faceOwner()[facei])
344 mesh().isInternalFace(facei)
345 && (nei !=
mesh().faceNeighbour()[facei])
347 || (newFace !=
mesh().faces()[facei])
352 Pout<<
"Modifying face " << facei
353 <<
" old vertices:" <<
mesh().
faces()[facei]
354 <<
" new vertices:" << newFace
355 <<
" new owner:" << own
356 <<
" new neighbour:" << nei
358 <<
" new zoneID:" <<
zoneID
359 <<
" new zoneFlip:" << zoneFlip
363 if ((nei == -1) || (own != -1 && own < nei))
387 newFace.reverseFace(),
403void Foam::meshCutAndRemove::copyFace
417 newFace[newFp++] =
f[fp];
419 fp = (fp + 1) %
f.
size();
421 newFace[newFp] =
f[fp];
429void Foam::meshCutAndRemove::splitFace
440 label startFp =
f.
find(v0);
445 <<
"Cannot find vertex (new numbering) " << v0
450 label endFp =
f.
find(v1);
455 <<
"Cannot find vertex (new numbering) " << v1
461 f0.setSize((endFp + 1 +
f.
size() - startFp) %
f.
size());
462 f1.setSize(
f.
size() - f0.size() + 2);
464 copyFace(
f, startFp, endFp, f0);
465 copyFace(
f, endFp, startFp, f1);
469Foam::face Foam::meshCutAndRemove::addEdgeCutsToFace(
const label facei)
const
473 face newFace(2 *
f.
size());
480 newFace[newFp++] =
f[fp];
485 EdgeMap<label>::const_iterator fnd =
486 addedPoints_.find(edge(
f[fp],
f[fp1]));
491 newFace[newFp++] = fnd.val();
495 newFace.setSize(newFp);
510 face newFace(2*loop.size());
516 label
cut = loop[fp];
520 label edgeI = getEdge(
cut);
524 label vertI = addedPoints_[
e];
526 newFace[newFacei++] = vertI;
531 label vertI = getVertex(
cut);
533 newFace[newFacei++] = vertI;
535 label nextCut = loop[loop.fcIndex(fp)];
537 if (!isEdge(nextCut))
540 label nextVertI = getVertex(nextCut);
547 EdgeMap<label>::const_iterator fnd =
548 addedPoints_.find(
mesh().edges()[edgeI]);
552 newFace[newFacei++] = fnd.val();
558 newFace.setSize(newFacei);
578 const label exposedPatchi,
586 addedFaces_.resize(cuts.
nLoops());
588 addedPoints_.clear();
589 addedPoints_.resize(cuts.
nLoops());
600 if (exposedPatchi < 0 || exposedPatchi >=
patches.
size())
603 <<
"Illegal exposed patch " << exposedPatchi
619 if (debug && findCutCell(cuts,
mesh().edgeCells()[edgeI]) == -1)
622 <<
"Problem: cut edge but none of the cells using it is\n"
623 <<
"edge:" << edgeI <<
" verts:" <<
e
628 label masterPointi =
e.start();
635 point newPt = weight*v1 + (1.0-weight)*v0;
650 addedPoints_.insert(
e, addedPointi);
654 Pout<<
"Added point " << addedPointi
656 << masterPointi <<
" of edge " << edgeI
657 <<
" vertices " <<
e <<
endl;
671 const labelList& loop = cellLoops[celli];
678 label
cut = loop[fp];
682 usedPoint[getVertex(
cut)] =
true;
686 const labelList& anchors = anchorPts[celli];
690 usedPoint[anchors[i]] =
true;
700 usedPoint[cPoints[i]] =
true;
711 const edge& fCut = iter.val();
719 label pointi = getVertex(
cut);
721 if (!usedPoint[pointi])
724 <<
"Problem: faceSplitCut not used by any loop"
725 <<
" or cell anchor point"
726 <<
"face:" << iter.key() <<
" point:" << pointi
738 if (!usedPoint[pointi])
741 <<
"Problem: point is marked as cut but"
742 <<
" not used by any loop"
743 <<
" or cell anchor point"
744 <<
"point:" << pointi
755 if (!usedPoint[pointi])
761 Pout<<
"Removing unused point " << pointi <<
endl;
774 const labelList& loop = cellLoops[celli];
778 if (cutPatch[celli] < 0 || cutPatch[celli] >=
patches.
size())
781 <<
"Illegal patch " << cutPatch[celli]
782 <<
" provided for cut cell " << celli
790 face newFace(loopToFace(celli, loop));
795 label masterPointi = findPatchFacePoint(newFace, exposedPatchi);
815 addedFaces_.insert(celli, addedFacei);
819 Pout<<
"Added splitting face " << newFace <<
" index:"
820 << addedFacei <<
" from masterPoint:" << masterPointi
821 <<
" to owner " << celli <<
" with anchors:"
838 writeCuts(
Pout, loop, weights);
859 const label facei = iter.key();
861 const edge& splitEdge = iter.val();
864 face newFace(addEdgeCutsToFace(facei));
867 label cut0 = splitEdge[0];
872 label edgeI = getEdge(cut0);
873 v0 = addedPoints_[
mesh().
edges()[edgeI]];
877 v0 = getVertex(cut0);
880 label cut1 = splitEdge[1];
884 label edgeI = getEdge(cut1);
885 v1 = addedPoints_[
mesh().
edges()[edgeI]];
889 v1 = getVertex(cut1);
894 splitFace(newFace, v0, v1, f0, f1);
900 if (
mesh().isInternalFace(facei))
908 <<
" own:" << own <<
" nei:" << nei
910 <<
" and f1:" << f1 <<
endl;
929 if (cellLoops[own].empty())
935 else if (isIn(splitEdge, cellLoops[own]))
939 if (firstCommon(f0, anchorPts[own]) != -1)
956 if (firstCommon(
f, anchorPts[own]) != -1)
976 if (cellLoops[nei].empty())
981 else if (isIn(splitEdge, cellLoops[nei]))
987 if (firstCommon(f0, anchorPts[nei]) != -1)
1003 if (firstCommon(
f, anchorPts[nei]) != -1)
1020 Pout<<
"f0 own:" << f0Own <<
" nei:" << f0Nei
1021 <<
" f1 own:" << f1Own <<
" nei:" << f1Nei
1039 bool modifiedFacei =
false;
1046 modFace(meshMod, facei, f0, f0Own, f0Nei,
patchID);
1047 modifiedFacei =
true;
1055 modFace(meshMod, facei, f0, f0Own, f0Nei,
patchID);
1056 modifiedFacei =
true;
1061 modFace(meshMod, facei, f0, f0Own, f0Nei, -1);
1062 modifiedFacei =
true;
1080 modFace(meshMod, facei, f1, f1Own, f1Nei,
patchID);
1081 modifiedFacei =
true;
1085 label masterPointi = findPatchFacePoint(f1,
patchID);
1107 modFace(meshMod, facei, f1, f1Own, f1Nei,
patchID);
1108 modifiedFacei =
true;
1112 label masterPointi = findPatchFacePoint(f1,
patchID);
1131 modFace(meshMod, facei, f1, f1Own, f1Nei, -1);
1132 modifiedFacei =
true;
1136 label masterPointi = findPatchFacePoint(f1, -1);
1138 addFace(meshMod, facei, masterPointi, f1, f1Own, f1Nei, -1);
1143 if (f0Own == -1 && f0Nei == -1 && !modifiedFacei)
1149 Pout<<
"Removed face " << facei <<
endl;
1153 faceUptodate[facei] =
true;
1166 if (edgeIsCut[edgeI])
1172 label facei = eFaces[i];
1174 if (!faceUptodate[facei])
1185 if (own == -1 && nei == -1)
1191 Pout<<
"Removed face " << facei <<
endl;
1197 face newFace(addEdgeCutsToFace(facei));
1201 Pout<<
"Added edge cuts to face " << facei
1203 <<
" newFace:" << newFace <<
endl;
1217 faceUptodate[facei] =
true;
1234 if (!faceUptodate[facei])
1240 if (own == -1 && nei == -1)
1246 Pout<<
"Removed face " << facei <<
endl;
1251 modFace(meshMod, facei, faces[facei], own, nei,
patchID);
1254 faceUptodate[facei] =
true;
1260 Pout<<
"meshCutAndRemove:" <<
nl
1261 <<
" cells split:" << cuts.
nLoops() <<
nl
1262 <<
" faces added:" << addedFaces_.size() <<
nl
1263 <<
" points added on edges:" << addedPoints_.size() <<
nl
1273 Map<label> newAddedFaces(addedFaces_.size());
1277 const label celli = iter.key();
1278 const label addedFacei = iter.val();
1283 if ((newCelli >= 0) && (newAddedFacei >= 0))
1288 && (newCelli != celli || newAddedFacei != addedFacei)
1291 Pout<<
"meshCutAndRemove::updateMesh :"
1292 <<
" updating addedFace for cell " << celli
1293 <<
" from " << addedFacei
1294 <<
" to " << newAddedFacei
1297 newAddedFaces.
insert(newCelli, newAddedFacei);
1302 addedFaces_.transfer(newAddedFaces);
1310 const edge&
e = iter.key();
1311 const label addedPointi = iter.val();
1317 if ((newStart >= 0) && (newEnd >= 0) && (newAddedPointi >= 0))
1319 edge newE(newStart, newEnd);
1324 && (
e != newE || newAddedPointi != addedPointi)
1327 Pout<<
"meshCutAndRemove::updateMesh :"
1328 <<
" updating addedPoints for edge " <<
e
1329 <<
" from " << addedPointi
1330 <<
" to " << newAddedPointi
1334 newAddedPoints.
insert(newE, newAddedPointi);
1339 addedPoints_.transfer(newAddedPoints);
Map from edge (expressed as its endpoints) to value. For easier forward declaration it is currently i...
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
A HashTable to objects of type <T> with a label key.
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
iterator end() noexcept
Return an iterator to end traversing the UList.
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
void size(const label n)
Older name for setAddressableSize.
label fcIndex(const label i) const noexcept
label size() const noexcept
The number of elements in the list.
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Description of cuts across cells.
const labelListList & cellLoops() const
For each cut cell the cut along the circumference.
const boolList & pointIsCut() const
Is mesh point cut.
const boolList & edgeIsCut() const
Is edge cut.
const labelListList & cellAnchorPoints() const
For each cut cell the points on the 'anchor' side of the cell.
label nLoops() const
Number of valid cell loops.
const scalarField & edgeWeight() const
If edge is cut gives weight (ratio between start() and end())
const Map< edge > & faceSplitCut() const
Gives for split face the two cuts that split the face into two.
Combines edge or vertex in single label. Used to specify cuts across cell circumference.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Smooth ATC in cells next to a set of patches supplied by type.
A face is a list of labels corresponding to mesh vertices.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
const labelList & reverseCellMap() const
Reverse cell map.
const labelList & reversePointMap() const
Reverse point map.
const labelList & reverseFaceMap() const
Reverse face map.
Like meshCutter but also removes non-anchor side of cell.
void setRefinement(const label exposedPatchi, const cellCuts &cuts, const labelList &cutPatch, polyTopoChange &meshMod)
Do actual cutting with cut description. Inserts mesh changes.
void updateMesh()
Update for new mesh topology.
A face addition data class. A face can be inflated either from a point or from another face and can e...
Class containing data for point addition.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Mesh consisting of general polyhedral cells.
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
virtual const faceList & faces() const
Return raw faces.
virtual const labelList & faceOwner() const
Return face owner.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
virtual const labelList & faceNeighbour() const
Return face neighbour.
virtual const pointField & points() const
Return raw points.
Class containing data for face removal.
Class containing data for point removal.
Direct mesh changes based on v1.3 polyTopoChange syntax.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const labelListList & cellPoints() const
const labelListList & pointFaces() const
const labelListList & edgeFaces() const
Patchify triangles based on orientation w.r.t other (triangulated or triangulatable) surfaces.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const polyBoundaryMesh & patches
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const labelIOList & zoneID
List< label > labelList
A List of labels.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
Ostream & endl(Ostream &os)
Add newline and flush stream.
List< labelList > labelListList
A List of labelList.
errorManip< error > abort(error &err)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
constexpr char nl
The newline '\n' character (0x0a)
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]
labelList pointLabels(nPoints, -1)
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.