51bool Foam::meshCutter::uses(
const labelList& elems1,
const labelList& elems2)
55 if (elems2.found(elems1[elemI]))
64bool Foam::meshCutter::isIn
70 label index = cuts.find(twoCuts[0]);
79 cuts[cuts.fcIndex(index)] == twoCuts[1]
80 || cuts[cuts.rcIndex(index)] == twoCuts[1]
87Foam::label Foam::meshCutter::findCutCell
95 label celli = cellLabels[labelI];
97 if (cuts.cellLoops()[celli].size())
106Foam::label Foam::meshCutter::findInternalFacePoint
119 label facei =
pFaces[pFacei];
121 if (
mesh().isInternalFace(facei))
140 const cellCuts& cuts,
153 if (cellLoops[own].size() && uses(
f, anchorPts[own]))
155 own = addedCells_[own];
160 if (
mesh().isInternalFace(facei))
164 if (cellLoops[nei].size() && uses(
f, anchorPts[nei]))
166 nei = addedCells_[nei];
172void Foam::meshCutter::getFaceInfo
182 if (!
mesh().isInternalFace(facei))
195 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
200void Foam::meshCutter::addFace
202 polyTopoChange& meshMod,
213 if ((nei == -1) || (own < nei))
218 Pout<<
"Adding face " << newFace
219 <<
" with new owner:" << own
220 <<
" with new neighbour:" << nei
223 <<
" zoneFlip:" << zoneFlip
249 Pout<<
"Adding (reversed) face " << newFace.reverseFace()
250 <<
" with new owner:" << nei
251 <<
" with new neighbour:" << own
254 <<
" zoneFlip:" << zoneFlip
262 newFace.reverseFace(),
278void Foam::meshCutter::modFace
280 polyTopoChange& meshMod,
293 (own !=
mesh().faceOwner()[facei])
295 mesh().isInternalFace(facei)
296 && (nei !=
mesh().faceNeighbour()[facei])
298 || (newFace !=
mesh().faces()[facei])
303 Pout<<
"Modifying face " << facei
304 <<
" old vertices:" <<
mesh().
faces()[facei]
305 <<
" new vertices:" << newFace
306 <<
" new owner:" << own
307 <<
" new neighbour:" << nei
308 <<
" new zoneID:" <<
zoneID
309 <<
" new zoneFlip:" << zoneFlip
313 if ((nei == -1) || (own < nei))
337 newFace.reverseFace(),
353void Foam::meshCutter::copyFace
367 newFace[newFp++] =
f[fp];
369 fp = (fp + 1) %
f.
size();
371 newFace[newFp] =
f[fp];
375void Foam::meshCutter::splitFace
386 label startFp =
f.
find(v0);
391 <<
"Cannot find vertex (new numbering) " << v0
396 label endFp =
f.
find(v1);
401 <<
"Cannot find vertex (new numbering) " << v1
407 f0.setSize((endFp + 1 +
f.
size() - startFp) %
f.
size());
408 f1.setSize(
f.
size() - f0.size() + 2);
410 copyFace(
f, startFp, endFp, f0);
411 copyFace(
f, endFp, startFp, f1);
415Foam::face Foam::meshCutter::addEdgeCutsToFace(
const label facei)
const
419 face newFace(2 *
f.
size());
426 newFace[newFp++] =
f[fp];
431 EdgeMap<label>::const_iterator fnd =
432 addedPoints_.find(edge(
f[fp],
f[fp1]));
437 newFace[newFp++] = fnd.val();
441 newFace.setSize(newFp);
453 face newFace(2*loop.size());
459 label
cut = loop[fp];
463 label edgeI = getEdge(
cut);
467 label vertI = addedPoints_[
e];
469 newFace[newFacei++] = vertI;
474 label vertI = getVertex(
cut);
476 newFace[newFacei++] = vertI;
478 label nextCut = loop[loop.fcIndex(fp)];
480 if (!isEdge(nextCut))
483 label nextVertI = getVertex(nextCut);
490 EdgeMap<label>::const_iterator fnd =
491 addedPoints_.find(
mesh().edges()[edgeI]);
495 newFace[newFacei++] = fnd.val();
501 newFace.setSize(newFacei);
529 addedCells_.resize(cuts.
nLoops());
532 addedFaces_.resize(cuts.
nLoops());
534 addedPoints_.clear();
535 addedPoints_.resize(cuts.
nLoops());
556 edgeOnCutCell[cEdges[i]] =
true;
564 if (cuts.
edgeIsCut()[edgeI] && !edgeOnCutCell[edgeI])
569 <<
"Problem: cut edge but none of the cells using"
571 <<
"edge:" << edgeI <<
" verts:" <<
e
590 label masterPointi =
e.start();
597 point newPt = weight*v1 + (1.0-weight)*v0;
612 addedPoints_.insert(
e, addedPointi);
616 Pout<<
"Added point " << addedPointi
618 << masterPointi <<
" of edge " << edgeI
619 <<
" vertices " <<
e <<
endl;
630 if (cellLoops[celli].size())
642 mesh().cellZones().whichZone(celli)
646 addedCells_.insert(celli, addedCelli);
650 Pout<<
"Added cell " << addedCells_[celli] <<
" to cell "
663 const labelList& loop = cellLoops[celli];
670 face newFace(loopToFace(celli, loop));
673 label masterPointi = findInternalFacePoint(anchorPts[celli]);
693 addedFaces_.insert(celli, addedFacei);
711 Pout<<
"Added splitting face " << newFace <<
" index:"
713 <<
" to owner " << celli
714 <<
" neighbour " << addedCells_[celli]
716 writeCuts(
Pout, loop, weights);
736 const label facei = iter.key();
738 const edge& splitEdge = iter.val();
741 face newFace(addEdgeCutsToFace(facei));
744 label cut0 = splitEdge[0];
749 label edgeI = getEdge(cut0);
750 v0 = addedPoints_[
mesh().
edges()[edgeI]];
754 v0 = getVertex(cut0);
757 label cut1 = splitEdge[1];
761 label edgeI = getEdge(cut1);
762 v1 = addedPoints_[
mesh().
edges()[edgeI]];
766 v1 = getVertex(cut1);
771 splitFace(newFace, v0, v1, f0, f1);
777 if (
mesh().isInternalFace(facei))
785 <<
" own:" << own <<
" nei:" << nei
787 <<
" and f1:" << f1 <<
endl;
803 if (cellLoops[own].empty())
808 else if (isIn(splitEdge, cellLoops[own]))
812 if (uses(f0, anchorPts[own]))
814 f0Owner = addedCells_[own];
820 f1Owner = addedCells_[own];
827 if (uses(
f, anchorPts[own]))
829 label newCelli = addedCells_[own];
841 label f0Neighbour = -1;
842 label f1Neighbour = -1;
846 if (cellLoops[nei].empty())
851 else if (isIn(splitEdge, cellLoops[nei]))
855 if (uses(f0, anchorPts[nei]))
857 f0Neighbour = addedCells_[nei];
863 f1Neighbour = addedCells_[nei];
870 if (uses(
f, anchorPts[nei]))
872 label newCelli = addedCells_[nei];
873 f0Neighbour = newCelli;
874 f1Neighbour = newCelli;
885 addFace(meshMod, facei, f0, f0Owner, f0Neighbour);
887 modFace(meshMod, facei, f1, f1Owner, f1Neighbour);
889 faceUptodate[facei] =
true;
902 if (edgeIsCut[edgeI])
908 label facei = eFaces[i];
910 if (!faceUptodate[facei])
913 face newFace(addEdgeCutsToFace(facei));
917 Pout<<
"Added edge cuts to face " << facei
919 <<
" newFace:" << newFace <<
endl;
926 modFace(meshMod, facei, newFace, own, nei);
928 faceUptodate[facei] =
true;
941 if (cellLoops[celli].size())
945 forAll(cllFaces, cllFacei)
947 label facei = cllFaces[cllFacei];
949 if (!faceUptodate[facei])
954 if (debug && (
f != addEdgeCutsToFace(facei)))
957 <<
"Problem: edges added to face which does "
958 <<
" not use a marked cut" <<
endl
959 <<
"facei:" << facei <<
endl
960 <<
"face:" <<
f <<
endl
961 <<
"newFace:" << addEdgeCutsToFace(facei)
978 faceUptodate[facei] =
true;
986 Pout<<
"meshCutter:" <<
nl
987 <<
" cells split:" << addedCells_.size() <<
nl
988 <<
" faces added:" << addedFaces_.size() <<
nl
989 <<
" points added on edges:" << addedPoints_.size() <<
nl
1002 Map<label> newAddedCells(addedCells_.size());
1006 const label celli = iter.key();
1007 const label addedCelli = iter.val();
1010 const label newAddedCelli = morphMap.
reverseCellMap()[addedCelli];
1012 if (newCelli >= 0 && newAddedCelli >= 0)
1017 && (newCelli != celli || newAddedCelli != addedCelli)
1020 Pout<<
"meshCutter::updateMesh :"
1021 <<
" updating addedCell for cell " << celli
1022 <<
" from " << addedCelli
1023 <<
" to " << newAddedCelli <<
endl;
1025 newAddedCells.
insert(newCelli, newAddedCelli);
1030 addedCells_.transfer(newAddedCells);
1034 Map<label> newAddedFaces(addedFaces_.size());
1038 const label celli = iter.key();
1039 const label addedFacei = iter.val();
1042 const label newAddedFacei = morphMap.
reverseFaceMap()[addedFacei];
1044 if ((newCelli >= 0) && (newAddedFacei >= 0))
1049 && (newCelli != celli || newAddedFacei != addedFacei)
1052 Pout<<
"meshCutter::updateMesh :"
1053 <<
" updating addedFace for cell " << celli
1054 <<
" from " << addedFacei
1055 <<
" to " << newAddedFacei
1058 newAddedFaces.
insert(newCelli, newAddedFacei);
1063 addedFaces_.transfer(newAddedFaces);
1071 const edge&
e = iter.key();
1072 const label addedPointi = iter.val();
1080 if ((newStart >= 0) && (newEnd >= 0) && (newAddedPointi >= 0))
1082 edge newE =
edge(newStart, newEnd);
1087 && (
e != newE || newAddedPointi != addedPointi)
1090 Pout<<
"meshCutter::updateMesh :"
1091 <<
" updating addedPoints for edge " <<
e
1092 <<
" from " << addedPointi
1093 <<
" to " << newAddedPointi
1097 newAddedPoints.
insert(newE, newAddedPointi);
1102 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 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 & 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.
void setRefinement(const cellCuts &cuts, polyTopoChange &meshMod)
Do actual cutting with cut description. Inserts mesh changes.
void updateMesh()
Update for new mesh topology.
Class containing data for cell addition.
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.
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.
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 & cellEdges() const
const labelListList & pointFaces() const
const labelListList & edgeFaces() const
const cellList & cells() const
Patchify triangles based on orientation w.r.t other (triangulated or triangulatable) surfaces.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const labelList nEdges(UPstream::listGatherValues< label >(aMesh.nEdges()))
const labelIOList & zoneID
#define WarningInFunction
Report a warning using Foam::Warning.
List< label > labelList
A List of labels.
Ostream & endl(Ostream &os)
Add newline and flush stream.
List< labelList > labelListList
A List of labelList.
errorManip< error > abort(error &err)
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
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.