52 Foam::label Foam::meshCutAndRemove::firstCommon
60 label index1 = elems2.find(elems1[elemI]);
72 bool Foam::meshCutAndRemove::isIn
78 label index = cuts.find(twoCuts[0]);
87 cuts[cuts.fcIndex(index)] == twoCuts[1]
88 || cuts[cuts.rcIndex(index)] == twoCuts[1]
95 Foam::label Foam::meshCutAndRemove::findCutCell
101 forAll(cellLabels, labelI)
103 label celli = cellLabels[labelI];
105 if (cuts.cellLoops()[celli].size())
114 Foam::label Foam::meshCutAndRemove::findInternalFacePoint
127 label facei =
pFaces[pFacei];
129 if (
mesh().isInternalFace(facei))
146 Foam::label Foam::meshCutAndRemove::findPatchFacePoint
149 const label exposedPatchi
157 label pointi =
f[fp];
176 void Foam::meshCutAndRemove::faceCells
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))
221 void Foam::meshCutAndRemove::getZoneInfo
236 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
241 void 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
267 <<
" zoneID:" << zoneID
268 <<
" zoneFlip:" << zoneFlip
294 Pout<<
"Adding (reversed) face " << newFace.reverseFace()
295 <<
" with new owner:" << nei
296 <<
" with new neighbour:" << own
298 <<
" anchor:" << masterPointi
299 <<
" zoneID:" << zoneID
300 <<
" zoneFlip:" << zoneFlip
308 newFace.reverseFace(),
325 void 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(),
403 void Foam::meshCutAndRemove::copyFace
417 newFace[newFp++] =
f[fp];
419 fp = (fp + 1) %
f.
size();
421 newFace[newFp] =
f[fp];
429 void 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);
469 Foam::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];
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);
566 Foam::meshCutAndRemove::meshCutAndRemove(
const polyMesh&
mesh)
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])
757 meshMod.
setAction(polyRemovePoint(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];
873 v0 = addedPoints_[
mesh().
edges()[edgeI]];
877 v0 = getVertex(cut0);
880 label cut1 = splitEdge[1];
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)
1145 meshMod.
setAction(polyRemoveFace(facei));
1149 Pout<<
"Removed face " << facei <<
endl;
1153 faceUptodate[facei] =
true;
1166 if (edgeIsCut[edgeI])
1172 label facei = eFaces[i];
1174 if (!faceUptodate[facei])
1182 faceCells(cuts, exposedPatchi, facei, own, nei,
patchID);
1185 if (own == -1 && nei == -1)
1187 meshMod.
setAction(polyRemoveFace(facei));
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])
1238 faceCells(cuts, exposedPatchi, facei, own, nei,
patchID);
1240 if (own == -1 && nei == -1)
1242 meshMod.
setAction(polyRemoveFace(facei));
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();
1280 const label newCelli = map.reverseCellMap()[celli];
1281 const label newAddedFacei = map.reverseFaceMap()[addedFacei];
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);
1306 EdgeMap<label> newAddedPoints(addedPoints_.size());
1310 const edge&
e = iter.key();
1311 const label addedPointi = iter.val();
1313 const label newStart = map.reversePointMap()[
e.start()];
1314 const label newEnd = map.reversePointMap()[
e.end()];
1315 const label newAddedPointi = map.reversePointMap()[addedPointi];
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);
void size(const label n)
Older name for setAddressableSize.
const boolList & edgeIsCut() const
Is edge cut.
A face is a list of labels corresponding to mesh vertices.
labelList pointLabels(nPoints, -1)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
virtual const labelList & faceNeighbour() const
Return face neighbour.
constexpr char nl
The newline '\n' character (0x0a)
bool empty() const noexcept
True if List is empty (ie, size() is zero)
Ostream & endl(Ostream &os)
Add newline and flush stream.
void updateMesh(const mapPolyMesh &)
Force recalculation of locally stored data on topological change.
List< labelList > labelListList
List of labelList.
Description of cuts across cells.
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
List< face > faceList
List of faces.
static label getEdge(List< DynamicList< label >> &pe, DynamicList< edge > &es, const label pointi, const label nextPointi)
const dimensionedScalar e
Elementary charge.
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
label size() const noexcept
The number of entries in the list.
virtual const labelList & faceOwner() const
Return face owner.
Combines edge or vertex in single label. Used to specify cuts across cell circumference.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
virtual const faceList & faces() const
Return raw faces.
label whichPatch(const label meshFacei) const
Return patch index for a given mesh face index. Uses binary search.
void setRefinement(const label exposedPatchi, const cellCuts &cuts, const labelList &cutPatch, polyTopoChange &meshMod)
Do actual cutting with cut description. Inserts mesh changes.
errorManip< error > abort(error &err)
int debug
Static debugging option.
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
defineTypeNameAndDebug(combustionModel, 0)
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]
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const scalarField & edgeWeight() const
If edge is cut gives weight (ratio between start() and end())
label nLoops() const
Number of valid cell loops.
const boolList & pointIsCut() const
Is mesh point cut.
const Map< edge > & faceSplitCut() const
Gives for split face the two cuts that split the face into two.
vector point
Point is a vector.
const labelListList & cellLoops() const
For each cut cell the cut along the circumference.
Direct mesh changes based on v1.3 polyTopoChange syntax.
const polyBoundaryMesh & patches
const labelListList & pointFaces() const
Mesh consisting of general polyhedral cells.
List< label > labelList
A List of labels.
List< bool > boolList
A List of bools.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
const labelListList & edgeFaces() const
const labelListList & cellPoints() const
forAllConstIters(mixture.phases(), phase)
const labelListList & cellAnchorPoints() const
For each cut cell the points on the 'anchor' side of the cell.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.