55 if (elems2.found(elems1[elemI]))
64 bool Foam::meshCutter::isIn
70 label index = cuts.find(twoCuts[0]);
79 cuts[cuts.fcIndex(index)] == twoCuts[1]
80 || cuts[cuts.rcIndex(index)] == twoCuts[1]
87 Foam::label Foam::meshCutter::findCutCell
95 label celli = cellLabels[labelI];
97 if (cuts.cellLoops()[celli].size())
106 Foam::label Foam::meshCutter::findInternalFacePoint
119 label facei =
pFaces[pFacei];
121 if (
mesh().isInternalFace(facei))
138 void Foam::meshCutter::faceCells
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];
172 void Foam::meshCutter::getFaceInfo
182 if (!
mesh().isInternalFace(facei))
195 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
200 void Foam::meshCutter::addFace
202 polyTopoChange& meshMod,
209 label
patchID, zoneID, zoneFlip;
211 getFaceInfo(facei,
patchID, zoneID, zoneFlip);
213 if ((nei == -1) || (own < nei))
218 Pout<<
"Adding face " << newFace
219 <<
" with new owner:" << own
220 <<
" with new neighbour:" << nei
222 <<
" zoneID:" << zoneID
223 <<
" zoneFlip:" << zoneFlip
249 Pout<<
"Adding (reversed) face " << newFace.reverseFace()
250 <<
" with new owner:" << nei
251 <<
" with new neighbour:" << own
253 <<
" zoneID:" << zoneID
254 <<
" zoneFlip:" << zoneFlip
262 newFace.reverseFace(),
278 void Foam::meshCutter::modFace
280 polyTopoChange& meshMod,
287 label
patchID, zoneID, zoneFlip;
289 getFaceInfo(facei,
patchID, zoneID, zoneFlip);
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(),
353 void Foam::meshCutter::copyFace
367 newFace[newFp++] =
f[fp];
369 fp = (fp + 1) %
f.
size();
371 newFace[newFp] =
f[fp];
375 void 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);
415 Foam::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];
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_.reserve(cuts.
nLoops());
532 addedFaces_.reserve(cuts.
nLoops());
534 addedPoints_.clear();
535 addedPoints_.reserve(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];
750 v0 = addedPoints_[
mesh().
edges()[edgeI]];
754 v0 = getVertex(cut0);
757 label cut1 = splitEdge[1];
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;
924 faceCells(cuts, facei, own, nei);
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)
967 faceCells(cuts, facei, own, nei);
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();
1009 const label newCelli = morphMap.reverseCellMap()[celli];
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();
1041 const label newCelli = morphMap.reverseCellMap()[celli];
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);
1067 EdgeMap<label> newAddedPoints(addedPoints_.size());
1071 const edge&
e = iter.key();
1072 const label addedPointi = iter.val();
1074 label newStart = morphMap.reversePointMap()[
e.start()];
1076 label newEnd = morphMap.reversePointMap()[
e.end()];
1078 label newAddedPointi = morphMap.reversePointMap()[addedPointi];
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);
const labelListList & cellEdges() const
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.
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.
const cellList & cells() const
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...
void updateMesh(const mapPolyMesh &)
Force recalculation of locally stored data on topological change.
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.
virtual const labelList & faceOwner() const
Return face owner.
Combines edge or vertex in single label. Used to specify cuts across cell circumference.
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.
errorManip< error > abort(error &err)
label find(const T &val) const
Find index of the first occurrence of the value.
void setRefinement(const cellCuts &cuts, polyTopoChange &meshMod)
Do actual cutting with cut description. Inserts mesh changes.
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 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.
#define WarningInFunction
Report a warning using Foam::Warning.
Direct mesh changes based on v1.3 polyTopoChange syntax.
const labelListList & pointFaces() const
Mesh consisting of general polyhedral cells.
List< label > labelList
A List of labels.
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
List< bool > boolList
A List of bools.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
const labelListList & edgeFaces() 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.