63Foam::label Foam::cellCuts::findPartIndex
65 const labelList& elems,
70 for (label i = 0; i < nElems; i++)
91 result[labels[labelI]] =
true;
108 result[labels[labelI]] = weights[labelI];
114Foam::label Foam::cellCuts::firstUnique
117 const Map<label>& map
122 if (!map.found(lst[i]))
133void Foam::cellCuts::syncProc()
173 const polyPatch& pp = pbm[patchi];
179 label facei = pp.start()+i;
182 const auto iter = faceSplitCut_.cfind(facei);
188 const edge& cuts = iter.val();
194 label edgei = getEdge(cuts[i]);
195 label index = fEdges.find(edgei);
196 relCuts[bFacei][i] = -index-1;
200 label index =
f.
find(getVertex(cuts[i]));
201 relCuts[bFacei][i] = index+1;
221 const polyPatch& pp = pbm[patchi];
227 label facei = pp.start()+i;
230 const edge& relCut = relCuts[bFacei];
231 if (relCut != edge(0, 0))
236 edge absoluteCut(0, 0);
241 label oppFp = -relCut[i]-1;
242 label fp =
f.
size()-1-oppFp;
243 absoluteCut[i] = edgeToEVert(fEdges[fp]);
247 label oppFp = relCut[i]-1;
254 absoluteCut[i] = vertToEVert(
f[fp]);
260 !faceSplitCut_.insert(facei, absoluteCut)
261 && faceSplitCut_[facei] != absoluteCut
265 <<
"Cut " << faceSplitCut_[facei]
267 <<
" of coupled patch " << pp.name()
268 <<
" is not consistent with coupled cut "
280void Foam::cellCuts::writeUncutOBJ
287 OFstream cutsStream(dir /
"cell_" +
name(celli) +
".obj");
290 <<
" to " << cutsStream.name() <<
nl;
302 OFstream cutStream(dir /
"cellCuts_" +
name(celli) +
".obj");
305 <<
" to " << cutStream.name() <<
nl;
311 label pointi = cPoints[i];
312 if (pointIsCut_[pointi])
324 label edgeI = cEdges[i];
326 if (edgeIsCut_[edgeI])
330 const scalar w = edgeWeight_[edgeI];
338void Foam::cellCuts::writeOBJ
347 OFstream cutsStream(dir /
"cell_" +
name(celli) +
".obj");
350 <<
" to " << cutsStream.name() <<
nl;
363 OFstream loopStream(dir /
"cellLoop_" +
name(celli) +
".obj");
366 <<
" to " << loopStream.name() <<
nl;
370 writeOBJ(loopStream, loopPoints, vertI);
374 OFstream anchorStream(dir /
"anchors_" +
name(celli) +
".obj");
377 <<
" to " << anchorStream.name() <<
endl;
386Foam::label Foam::cellCuts::edgeEdgeToFace
397 label facei = cFaces[cFacei];
401 if (fEdges.found(edgeA) && fEdges.found(edgeB))
411 <<
"cellCuts : Cannot find face on cell "
412 << celli <<
" that has both edges " << edgeA <<
' ' << edgeB <<
endl
413 <<
"faces : " << cFaces <<
endl
416 <<
"Marking the loop across this cell as invalid" <<
endl;
422Foam::label Foam::cellCuts::edgeVertexToFace
433 label facei = cFaces[cFacei];
439 if (fEdges.found(edgeI) &&
f.
found(vertI))
446 <<
"cellCuts : Cannot find face on cell "
447 << celli <<
" that has both edge " << edgeI <<
" and vertex "
449 <<
"faces : " << cFaces <<
endl
451 <<
"Marking the loop across this cell as invalid" <<
endl;
457Foam::label Foam::cellCuts::vertexVertexToFace
468 label facei = cFaces[cFacei];
479 <<
"cellCuts : Cannot find face on cell "
480 << celli <<
" that has vertex " << vertA <<
" and vertex "
482 <<
"faces : " << cFaces <<
endl
483 <<
"Marking the loop across this cell as invalid" <<
endl;
489void Foam::cellCuts::calcFaceCuts()
const
500 auto& faceCuts = *faceCutsPtr_;
502 for (label facei = 0; facei <
mesh().
nFaces(); facei++)
504 const face&
f = faces[facei];
533 && !edgeIsCut_[
findEdge(facei, v0, vMin1)]
536 cuts[cutI++] = vertToEVert(v0);
553 label edgeI =
findEdge(facei, v0, v1);
555 if (edgeIsCut_[edgeI] && !pointIsCut_[v0])
557 cuts[cutI++] = edgeToEVert(edgeI);
574 bool allVerticesCut =
true;
584 label edgeI =
findEdge(facei, v0, v1);
588 cuts[cutI++] = vertToEVert(v0);
593 allVerticesCut =
false;
596 if (edgeIsCut_[edgeI])
598 cuts[cutI++] = edgeToEVert(edgeI);
607 if (verbose_ || debug)
610 <<
"Face " << facei <<
" vertices " <<
f
611 <<
" has all its vertices cut. Not cutting face." <<
endl;
618 if (cutI > 1 && cuts[cutI-1] == cuts[0])
627Foam::label Foam::cellCuts::findEdge
640 const edge&
e = edges[fEdges[i]];
644 (
e[0] == v0 &&
e[1] == v1)
645 || (
e[0] == v1 &&
e[1] == v0)
656Foam::label Foam::cellCuts::loopFace
662 const cell& cFaces =
mesh().
cells()[celli];
666 label facei = cFaces[cFacei];
671 bool allOnFace =
true;
679 if (!fEdges.found(getEdge(
cut)))
707bool Foam::cellCuts::walkPoint
710 const label startCut,
712 const label exclude0,
713 const label exclude1,
715 const label otherCut,
721 label vertI = getVertex(otherCut);
727 label otherFacei =
pFaces[pFacei];
731 otherFacei != exclude0
732 && otherFacei != exclude1
736 label oldNVisited = nVisited;
755 nVisited = oldNVisited;
762bool Foam::cellCuts::crossEdge
765 const label startCut,
767 const label otherCut,
774 label edgeI = getEdge(otherCut);
779 label oldNVisited = nVisited;
800 nVisited = oldNVisited;
807bool Foam::cellCuts::addCut
816 if (findPartIndex(visited, nVisited,
cut) != -1)
820 truncVisited.setSize(nVisited);
822 if (verbose_ || debug)
824 Pout<<
"For cell " << celli <<
" : trying to add duplicate cut "
827 writeCuts(
Pout, cuts, loopWeights(cuts));
830 writeCuts(
Pout, truncVisited, loopWeights(truncVisited));
838 visited[nVisited++] =
cut;
845bool Foam::cellCuts::walkFace
848 const label startCut,
853 label& beforeLastCut,
858 const labelList& fCuts = faceCuts()[facei];
860 if (fCuts.size() < 2)
866 if (fCuts.size() == 2)
870 if (!addCut(celli,
cut, nVisited, visited))
882 if (!addCut(celli,
cut, nVisited, visited))
899 for (label i = 0; i < fCuts.size()-1; i++)
901 if (!addCut(celli, fCuts[i], nVisited, visited))
906 beforeLastCut = fCuts[fCuts.size()-2];
907 lastCut = fCuts[fCuts.size()-1];
909 else if (fCuts[fCuts.size()-1] ==
cut)
911 for (label i = fCuts.size()-1; i >= 1; --i)
913 if (!addCut(celli, fCuts[i], nVisited, visited))
918 beforeLastCut = fCuts[1];
923 if (verbose_ || debug)
926 <<
"In middle of cut. cell:" << celli <<
" face:" << facei
927 <<
" cuts:" << fCuts <<
" current cut:" <<
cut <<
endl;
938bool Foam::cellCuts::walkCell
941 const label startCut,
951 label beforeLastCut = -1;
956 Pout<<
"For cell:" << celli <<
" walked across face " << facei
959 writeCuts(
Pout, cuts, loopWeights(cuts));
983 Pout<<
" to last cut ";
985 writeCuts(
Pout, cuts, loopWeights(cuts));
990 if (lastCut == startCut)
998 truncVisited.setSize(nVisited);
1000 Pout<<
"For cell " << celli <<
" : found closed path:";
1001 writeCuts(
Pout, truncVisited, loopWeights(truncVisited));
1002 Pout<<
" closed by " << lastCut <<
endl;
1030 if (isEdge(beforeLastCut))
1032 if (isEdge(lastCut))
1066 if (isEdge(lastCut))
1087 getVertex(beforeLastCut),
1131void Foam::cellCuts::calcCellLoops(
const labelList& cutCells)
1146 forAll(allFaceCuts, facei)
1148 const labelList& fCuts = allFaceCuts[facei];
1150 if (fCuts.size() ==
mesh().faces()[facei].size())
1155 if (
mesh().isInternalFace(facei))
1160 else if (fCuts.size() >= 2)
1165 if (
mesh().isInternalFace(facei))
1179 label celli = cutCells[i];
1181 bool validLoop =
false;
1184 if (nCutFaces[celli] >= 1)
1190 Pout<<
"cell:" << celli <<
" cut faces:" <<
endl;
1193 label facei = cFaces[i];
1194 const labelList& fCuts = allFaceCuts[facei];
1196 Pout<<
" face:" << facei <<
" cuts:";
1197 writeCuts(
Pout, fCuts, loopWeights(fCuts));
1207 label facei = cFaces[cFacei];
1209 const labelList& fCuts = allFaceCuts[facei];
1215 if (fCuts.size() >= 2)
1222 Pout<<
"cell:" << celli
1223 <<
" start walk at face:" << facei
1226 writeCuts(
Pout, cuts, loopWeights(cuts));
1261 for (label i = 0; i < nVisited; i++)
1263 loop[i] = visited[i];
1270 if (verbose_ || debug)
1272 Pout<<
"calcCellLoops(const labelList&) :"
1273 <<
" did not find valid"
1274 <<
" loop for cell " << celli <<
endl;
1276 writeUncutOBJ(
".", celli);
1278 cellLoops_[celli].clear();
1286 cellLoops_[celli].clear();
1292void Foam::cellCuts::walkEdges
1298 Map<label>& edgeStatus,
1299 Map<label>& pointStatus
1302 if (pointStatus.insert(pointi, status))
1310 label edgeI = pEdges[pEdgeI];
1315 && edgeStatus.insert(edgeI, status)
1320 label v2 =
mesh().
edges()[edgeI].otherVertex(pointi);
1322 walkEdges(celli, v2, status, edgeStatus, pointStatus);
1341 label pointi = cellPoints[i];
1345 !anchorPoints.
found(pointi)
1346 && !loop.
found(vertToEVert(pointi))
1349 newElems[newElemI++] = pointi;
1359bool Foam::cellCuts::loopAnchorConsistent
1369 const vector areaNorm =
f.areaNormal(loopPts);
1370 const point ctr =
f.centre(loopPts);
1376 for (
const label pointi : anchorPoints)
1380 avg /= anchorPoints.
size();
1383 if (((avg - ctr) & areaNorm) > 0)
1392bool Foam::cellCuts::calcAnchors
1405 const cell& cFaces =
mesh().
cells()[celli];
1413 Map<label> pointStatus(2*cPoints.size());
1414 Map<label> edgeStatus(2*cEdges.size());
1419 label
cut = loop[i];
1423 edgeStatus.insert(getEdge(
cut), 0);
1427 pointStatus.insert(getVertex(
cut), 0);
1434 label edgeI = cEdges[i];
1435 const edge&
e = edges[edgeI];
1437 if (pointStatus.found(
e[0]) && pointStatus.found(
e[1]))
1439 edgeStatus.insert(edgeI, 0);
1445 label uncutIndex = firstUnique(cPoints, pointStatus);
1447 if (uncutIndex == -1)
1449 if (verbose_ || debug)
1452 <<
"Invalid loop " << loop <<
" for cell " << celli <<
endl
1453 <<
"Can not find point on cell which is not cut by loop."
1463 walkEdges(celli, cPoints[uncutIndex], 1, edgeStatus, pointStatus);
1466 uncutIndex = firstUnique(cPoints, pointStatus);
1468 if (uncutIndex == -1)
1472 if (verbose_ || debug)
1475 <<
"Invalid loop " << loop <<
" for cell " << celli <<
endl
1476 <<
"All vertices of cell are either in loop or in anchor set"
1487 walkEdges(celli, cPoints[uncutIndex], 2, edgeStatus, pointStatus);
1490 DynamicList<label> connectedPoints(cPoints.size());
1491 DynamicList<label> otherPoints(cPoints.size());
1495 if (iter.val() == 1)
1497 connectedPoints.append(iter.key());
1499 else if (iter.val() == 2)
1501 otherPoints.append(iter.key());
1504 connectedPoints.shrink();
1505 otherPoints.shrink();
1508 uncutIndex = firstUnique(cPoints, pointStatus);
1510 if (uncutIndex != -1)
1512 if (verbose_ || debug)
1515 <<
"Invalid loop " << loop <<
" for cell " << celli
1516 <<
" since it splits the cell into more than two cells" <<
endl;
1518 writeOBJ(
".", celli, loopPts, connectedPoints);
1531 const label pointi = iter.key();
1535 if (iter.val() == 1)
1541 connectedFaces.insert(
pFaces[pFacei]);
1545 else if (iter.val() == 2)
1551 otherFaces.insert(
pFaces[pFacei]);
1557 if (connectedFaces.size() < 3)
1559 if (verbose_ || debug)
1562 <<
"Invalid loop " << loop <<
" for cell " << celli
1563 <<
" since would have too few faces on one side." <<
nl
1564 <<
"All faces:" << cFaces <<
endl;
1566 writeOBJ(
".", celli, loopPts, connectedPoints);
1572 if (otherFaces.size() < 3)
1574 if (verbose_ || debug)
1577 <<
"Invalid loop " << loop <<
" for cell " << celli
1578 <<
" since would have too few faces on one side." <<
nl
1579 <<
"All faces:" << cFaces <<
endl;
1581 writeOBJ(
".", celli, loopPts, otherPoints);
1595 label facei = cFaces[i];
1599 bool hasSet1 =
false;
1600 bool hasSet2 =
false;
1602 label prevStat = pointStatus[
f[0]];
1607 label pStat = pointStatus[v0];
1609 if (pStat == prevStat)
1612 else if (pStat == 0)
1616 else if (pStat == 1)
1621 if (verbose_ || debug)
1624 <<
"Invalid loop " << loop <<
" for cell "
1626 <<
" since face " <<
f <<
" would be split into"
1627 <<
" more than two faces" <<
endl;
1629 writeOBJ(
".", celli, loopPts, otherPoints);
1636 else if (pStat == 2)
1641 if (verbose_ || debug)
1644 <<
"Invalid loop " << loop <<
" for cell "
1646 <<
" since face " <<
f <<
" would be split into"
1647 <<
" more than two faces" <<
endl;
1649 writeOBJ(
".", celli, loopPts, otherPoints);
1665 label v1 =
f.nextLabel(fp);
1666 label edgeI =
findEdge(facei, v0, v1);
1668 label eStat = edgeStatus[edgeI];
1670 if (eStat == prevStat)
1673 else if (eStat == 0)
1677 else if (eStat == 1)
1682 if (verbose_ || debug)
1685 <<
"Invalid loop " << loop <<
" for cell "
1687 <<
" since face " <<
f <<
" would be split into"
1688 <<
" more than two faces" <<
endl;
1690 writeOBJ(
".", celli, loopPts, otherPoints);
1698 else if (eStat == 2)
1703 if (verbose_ || debug)
1706 <<
"Invalid loop " << loop <<
" for cell "
1708 <<
" since face " <<
f <<
" would be split into"
1709 <<
" more than two faces" <<
endl;
1711 writeOBJ(
".", celli, loopPts, otherPoints);
1727 bool loopOk = loopAnchorConsistent(celli, loopPts, connectedPoints);
1732 bool otherLoopOk = loopAnchorConsistent(celli, loopPts, otherPoints);
1734 if (loopOk == otherLoopOk)
1740 <<
" For cell:" << celli
1741 <<
" achorpoints and nonanchorpoints are geometrically"
1742 <<
" on same side!" <<
endl
1743 <<
"cellPoints:" << cPoints <<
endl
1744 <<
"loop:" << loop <<
endl
1745 <<
"anchors:" << connectedPoints <<
endl
1746 <<
"otherPoints:" << otherPoints <<
endl;
1748 writeOBJ(
".", celli, loopPts, connectedPoints);
1755 anchorPoints.transfer(connectedPoints);
1756 connectedPoints.clear();
1760 anchorPoints.transfer(otherPoints);
1761 otherPoints.clear();
1777 loopPts[fp] = coord(loop[fp], loopWeights[fp]);
1789 label
cut = loop[fp];
1793 label edgeI = getEdge(
cut);
1795 weights[fp] = edgeWeight_[edgeI];
1799 weights[fp] = -GREAT;
1806bool Foam::cellCuts::validEdgeLoop
1814 label
cut = loop[fp];
1818 label edgeI = getEdge(
cut);
1821 if (edgeIsCut_[edgeI])
1827 mag(loopWeights[fp] - edgeWeight_[edgeI])
1841Foam::label Foam::cellCuts::countFaceCuts
1858 label vertI =
f[fp];
1864 || loop.found(vertToEVert(vertI))
1876 label edgeI = fEdges[fEdgeI];
1882 || loop.found(edgeToEVert(edgeI))
1893bool Foam::cellCuts::conservativeValidLoop
1900 if (loop.size() < 2)
1907 if (isEdge(loop[cutI]))
1909 label edgeI = getEdge(loop[cutI]);
1911 if (edgeIsCut_[edgeI])
1920 if (pointIsCut_[
e.start()] || pointIsCut_[
e.
end()])
1931 label nCuts = countFaceCuts(eFaces[eFacei], loop);
1944 label vertI = getVertex(loop[cutI]);
1946 if (!pointIsCut_[vertI])
1955 label edgeI = pEdges[pEdgeI];
1957 if (edgeIsCut_[edgeI])
1968 label nCuts = countFaceCuts(
pFaces[pFacei], loop);
1984bool Foam::cellCuts::validLoop
1990 Map<edge>& newFaceSplitCut,
1999 if (loop.size() < 2)
2008 if (!conservativeValidLoop(celli, loop))
2010 Info <<
"Invalid conservative loop: " << loop <<
endl;
2017 label
cut = loop[fp];
2018 label nextCut = loop[(fp+1) % loop.size()];
2021 label meshFacei = -1;
2025 label edgeI = getEdge(
cut);
2029 if (isEdge(nextCut))
2032 label nextEdgeI = getEdge(nextCut);
2035 meshFacei = edgeEdgeToFace(celli, edgeI, nextEdgeI);
2037 if (meshFacei == -1)
2047 label nextVertI = getVertex(nextCut);
2050 if (
e.start() != nextVertI &&
e.
end() != nextVertI)
2053 meshFacei = edgeVertexToFace(celli, edgeI, nextVertI);
2055 if (meshFacei == -1)
2066 label vertI = getVertex(
cut);
2068 if (isEdge(nextCut))
2071 label nextEdgeI = getEdge(nextCut);
2073 const edge& nextE =
mesh().
edges()[nextEdgeI];
2075 if (nextE.start() != vertI && nextE.end() != vertI)
2078 meshFacei = edgeVertexToFace(celli, nextEdgeI, vertI);
2080 if (meshFacei == -1)
2089 label nextVertI = getVertex(nextCut);
2094 meshFacei = vertexVertexToFace(celli, vertI, nextVertI);
2096 if (meshFacei == -1)
2104 if (meshFacei != -1)
2108 edge cutEdge(
cut, nextCut);
2110 const auto iter = faceSplitCut_.cfind(meshFacei);
2115 newFaceSplitCut.insert(meshFacei, cutEdge);
2120 if (iter.val() != cutEdge)
2129 label faceContainingLoop = loopFace(celli, loop);
2131 if (faceContainingLoop != -1)
2133 if (verbose_ || debug)
2136 <<
"Found loop on cell " << celli <<
" with all points"
2137 <<
" on face " << faceContainingLoop <<
endl;
2151 loopPoints(loop, loopWeights),
2157void Foam::cellCuts::setFromCellLoops()
2160 pointIsCut_ =
false;
2164 faceSplitCut_.clear();
2166 forAll(cellLoops_, celli)
2168 const labelList& loop = cellLoops_[celli];
2173 Map<edge> faceSplitCuts(loop.size());
2190 if (verbose_ || debug)
2193 <<
"Illegal loop " << loop
2194 <<
" when recreating cut-addressing"
2195 <<
" from existing cellLoops for cell " << celli
2199 cellLoops_[celli].clear();
2200 cellAnchorPoints_[celli].clear();
2205 cellAnchorPoints_[celli].transfer(anchorPoints);
2210 faceSplitCut_.insert(iter.key(), iter.val());
2216 label
cut = loop[cutI];
2220 edgeIsCut_[getEdge(
cut)] =
true;
2224 pointIsCut_[getVertex(
cut)] =
true;
2232 forAll(edgeIsCut_, edgeI)
2234 if (!edgeIsCut_[edgeI])
2237 edgeWeight_[edgeI] = -GREAT;
2243bool Foam::cellCuts::setFromCellLoop
2256 OFstream str(
"last_cell.obj");
2258 str<<
"# edges of cell " << celli <<
nl;
2270 OFstream loopStr(
"last_loop.obj");
2272 loopStr<<
"# looppoints for cell " << celli <<
nl;
2274 pointField pointsOfLoop = loopPoints(loop, loopWeights);
2285 str <<
' ' << i + 1;
2287 str <<
' ' << 1 <<
nl;
2290 bool okLoop =
false;
2292 if (validEdgeLoop(loop, loopWeights))
2295 Map<edge> faceSplitCuts(loop.size());
2300 validLoop(celli, loop, loopWeights, faceSplitCuts, anchorPoints);
2305 cellLoops_[celli] = loop;
2306 cellAnchorPoints_[celli].transfer(anchorPoints);
2311 faceSplitCut_.insert(iter.key(), iter.val());
2318 const label
cut = loop[cutI];
2322 const label edgeI = getEdge(
cut);
2324 edgeIsCut_[edgeI] =
true;
2326 edgeWeight_[edgeI] = loopWeights[cutI];
2330 const label vertI = getVertex(
cut);
2332 pointIsCut_[vertI] =
true;
2342void Foam::cellCuts::setFromCellLoops
2346 const List<scalarField>& cellLoopWeights
2350 pointIsCut_ =
false;
2354 forAll(cellLabels, cellLabelI)
2356 const label celli = cellLabels[cellLabelI];
2358 const labelList& loop = cellLoops[cellLabelI];
2362 const scalarField& loopWeights = cellLoopWeights[cellLabelI];
2364 if (setFromCellLoop(celli, loop, loopWeights))
2371 cellLoops_[celli].
clear();
2378void Foam::cellCuts::setFromCellCutter
2380 const cellLooper& cellCutter,
2381 const List<refineCell>& refCells
2385 pointIsCut_ =
false;
2394 DynamicList<label> invalidCutCells(2);
2395 DynamicList<labelList> invalidCutLoops(2);
2396 DynamicList<scalarField> invalidCutLoopWeights(2);
2399 forAll(refCells, refCelli)
2401 const refineCell& refCell = refCells[refCelli];
2403 const label celli = refCell.cellNo();
2405 const vector& refDir = refCell.direction();
2427 if (setFromCellLoop(celli, cellLoop, cellLoopWeights))
2433 cellLoops_[celli].clear();
2436 <<
"Found loop on cell " << celli
2437 <<
" that resulted in an unexpected bad cut." <<
nl
2438 <<
" Suggestions:" <<
nl
2439 <<
" - Turn on the debug switch for 'cellCuts' to get"
2440 <<
" geometry files that identify this cell." <<
nl
2441 <<
" - Also keep in mind to check the defined"
2442 <<
" reference directions, as these are most likely the"
2443 <<
" origin of the problem."
2449 invalidCutCells.append(celli);
2450 invalidCutLoops.append(cellLoop);
2451 invalidCutLoopWeights.append(cellLoopWeights);
2458 cellLoops_[celli].clear();
2462 if (debug && invalidCutCells.size())
2464 invalidCutCells.shrink();
2465 invalidCutLoops.shrink();
2466 invalidCutLoopWeights.shrink();
2468 fileName cutsFile(
"invalidLoopCells.obj");
2470 Pout<<
"cellCuts : writing inValidLoops cells to " << cutsFile <<
endl;
2472 OFstream cutsStream(cutsFile);
2483 fileName loopsFile(
"invalidLoops.obj");
2485 Pout<<
"cellCuts : writing inValidLoops loops to " << loopsFile <<
endl;
2487 OFstream loopsStream(loopsFile);
2491 forAll(invalidCutLoops, i)
2496 loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
2504void Foam::cellCuts::setFromCellCutter
2506 const cellLooper& cellCutter,
2508 const List<plane>& cellCutPlanes
2512 pointIsCut_ =
false;
2521 DynamicList<label> invalidCutCells(2);
2522 DynamicList<labelList> invalidCutLoops(2);
2523 DynamicList<scalarField> invalidCutLoopWeights(2);
2528 const label celli = cellLabels[i];
2549 if (setFromCellLoop(celli, cellLoop, cellLoopWeights))
2555 cellLoops_[celli].clear();
2560 invalidCutCells.append(celli);
2561 invalidCutLoops.append(cellLoop);
2562 invalidCutLoopWeights.append(cellLoopWeights);
2569 cellLoops_[celli].clear();
2573 if (debug && invalidCutCells.size())
2575 invalidCutCells.shrink();
2576 invalidCutLoops.shrink();
2577 invalidCutLoopWeights.shrink();
2579 fileName cutsFile(
"invalidLoopCells.obj");
2581 Pout<<
"cellCuts : writing inValidLoops cells to " << cutsFile <<
endl;
2583 OFstream cutsStream(cutsFile);
2594 fileName loopsFile(
"invalidLoops.obj");
2596 Pout<<
"cellCuts : writing inValidLoops loops to " << loopsFile <<
endl;
2598 OFstream loopsStream(loopsFile);
2602 forAll(invalidCutLoops, i)
2607 loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
2615void Foam::cellCuts::orientPlanesAndLoops()
2618 forAll(cellLoops_, celli)
2620 const labelList& loop = cellLoops_[celli];
2622 if (loop.size() && cellAnchorPoints_[celli].empty())
2630 cellAnchorPoints_[celli]
2637 Pout<<
"cellAnchorPoints:" <<
endl;
2639 forAll(cellAnchorPoints_, celli)
2641 if (cellLoops_[celli].size())
2643 if (cellAnchorPoints_[celli].empty())
2646 <<
"No anchor points for cut cell " << celli <<
endl
2652 Pout<<
" cell:" << celli <<
" anchored at "
2653 << cellAnchorPoints_[celli] <<
endl;
2661 forAll(cellLoops_, celli)
2663 if (cellLoops_[celli].size())
2671void Foam::cellCuts::calcLoopsAndAddressing(
const labelList& cutCells)
2674 forAll(edgeIsCut_, edgeI)
2676 if (edgeIsCut_[edgeI])
2678 scalar weight = edgeWeight_[edgeI];
2680 if (weight < 0 || weight > 1)
2683 <<
"Weight out of range [0,1]. Edge " << edgeI
2691 edgeWeight_[edgeI] = -GREAT;
2697 calcCellLoops(cutCells);
2702 forAll(cellLoops_, celli)
2704 const labelList& loop = cellLoops_[celli];
2708 Pout<<
"cell:" << celli <<
" ";
2709 writeCuts(
Pout, loop, loopWeights(loop));
2721void Foam::cellCuts::check()
const
2724 forAll(edgeIsCut_, edgeI)
2726 if (edgeIsCut_[edgeI])
2736 <<
"edge:" << edgeI <<
" vertices:"
2738 <<
" weight:" << edgeWeight_[edgeI] <<
" should have been"
2739 <<
" snapped to one of its endpoints"
2745 if (edgeWeight_[edgeI] > - 1)
2748 <<
"edge:" << edgeI <<
" vertices:"
2750 <<
" weight:" << edgeWeight_[edgeI] <<
" is not cut but"
2751 <<
" its weight is not marked invalid"
2758 forAll(cellLoops_, celli)
2760 const labelList& loop = cellLoops_[celli];
2764 label
cut = loop[i];
2768 (isEdge(
cut) && !edgeIsCut_[getEdge(
cut)])
2769 || (!isEdge(
cut) && !pointIsCut_[getVertex(
cut)])
2773 writeCuts(
Pout, cuts, loopWeights(cuts));
2776 <<
"cell:" << celli <<
" loop:"
2778 <<
" cut:" <<
cut <<
" is not marked as cut"
2785 forAll(cellLoops_, celli)
2787 const labelList& anchors = cellAnchorPoints_[celli];
2789 const labelList& loop = cellLoops_[celli];
2791 if (loop.size() && anchors.empty())
2794 <<
"cell:" << celli <<
" loop:" << loop
2795 <<
" has no anchor points"
2802 label
cut = loop[i];
2807 && anchors.found(getVertex(
cut))
2811 <<
"cell:" << celli <<
" loop:" << loop
2812 <<
" anchor points:" << anchors
2813 <<
" anchor:" << getVertex(
cut) <<
" is part of loop"
2824 forAll(cellLoops_, celli)
2826 cellIsCut[celli] = cellLoops_[celli].size();
2833 const label facei = iter.key();
2835 if (
mesh().isInternalFace(facei))
2840 if (cellLoops_[own].empty() && cellLoops_[nei].empty())
2843 <<
"Internal face:" << facei <<
" cut by " << iter.val()
2844 <<
" has owner:" << own
2845 <<
" and neighbour:" << nei
2846 <<
" that are both uncut"
2856 if (cellLoops_[own].empty() && !nbrCellIsCut[bFacei])
2859 <<
"Boundary face:" << facei <<
" cut by " << iter.val()
2860 <<
" has owner:" << own
2884 edgeIsCut_(expand(
mesh.
nEdges(), meshEdges)),
2885 edgeWeight_(expand(
mesh.
nEdges(), meshEdges, meshEdgeWeights)),
2886 faceSplitCut_(cutCells.size()),
2887 cellLoops_(
mesh.nCells()),
2889 cellAnchorPoints_(
mesh.nCells())
2893 Pout<<
"cellCuts : constructor from cut verts and edges" <<
endl;
2896 calcLoopsAndAddressing(cutCells);
2899 orientPlanesAndLoops();
2910 Pout<<
"cellCuts : leaving constructor from cut verts and edges"
2928 edgeIsCut_(expand(
mesh.
nEdges(), meshEdges)),
2929 edgeWeight_(expand(
mesh.
nEdges(), meshEdges, meshEdgeWeights)),
2931 cellLoops_(
mesh.nCells()),
2933 cellAnchorPoints_(
mesh.nCells())
2940 Pout<<
"cellCuts : constructor from cellLoops" <<
endl;
2949 orientPlanesAndLoops();
2960 Pout<<
"cellCuts : leaving constructor from cellLoops" <<
endl;
2979 faceSplitCut_(cellLabels.size()),
2980 cellLoops_(
mesh.nCells()),
2982 cellAnchorPoints_(
mesh.nCells())
2986 Pout<<
"cellCuts : constructor from cellLoops" <<
endl;
2991 setFromCellLoops(cellLabels,
cellLoops, cellEdgeWeights);
2997 orientPlanesAndLoops();
3008 Pout<<
"cellCuts : leaving constructor from cellLoops" <<
endl;
3026 faceSplitCut_(refCells.size()),
3027 cellLoops_(
mesh.nCells()),
3029 cellAnchorPoints_(
mesh.nCells())
3033 Pout<<
"cellCuts : constructor from cellCutter" <<
endl;
3038 setFromCellCutter(cellCutter, refCells);
3044 orientPlanesAndLoops();
3055 Pout<<
"cellCuts : leaving constructor from cellCutter" <<
endl;
3074 faceSplitCut_(cellLabels.size()),
3075 cellLoops_(
mesh.nCells()),
3077 cellAnchorPoints_(
mesh.nCells())
3081 Pout<<
"cellCuts : constructor from cellCutter with prescribed plane"
3087 setFromCellCutter(cellCutter, cellLabels, cutPlanes);
3093 orientPlanesAndLoops();
3104 Pout<<
"cellCuts : leaving constructor from cellCutter with prescribed"
3105 <<
" plane" <<
endl;
3125 pointIsCut_(pointIsCut),
3126 edgeIsCut_(edgeIsCut),
3127 edgeWeight_(edgeWeight),
3128 faceSplitCut_(faceSplitCut),
3129 cellLoops_(cellLoops),
3131 cellAnchorPoints_(cellAnchorPoints)
3135 Pout<<
"cellCuts : constructor from components" <<
endl;
3136 Pout<<
"cellCuts : leaving constructor from components" <<
endl;
3145 faceCutsPtr_.reset(
nullptr);
3151 const labelList& loop = cellLoops_[celli];
3157 const label
cut = loop[fp];
3161 label edgeI = getEdge(
cut);
3163 loopPts[fp] = coord(
cut, edgeWeight_[edgeI]);
3167 loopPts[fp] = coord(
cut, -GREAT);
3181 cellAnchorPoints_[celli] =
3184 mesh().cellPoints()[celli],
3185 cellAnchorPoints_[celli],
3199void Foam::cellCuts::writeOBJ
3206 label startVertI = vertI;
3210 const point& pt = loopPts[fp];
3212 os <<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z() <<
endl;
3220 os <<
' ' << startVertI + fp + 1;
3230 forAll(cellLoops_, celli)
3232 if (cellLoops_[celli].size())
3234 writeOBJ(
os, loopPoints(celli), vertI);
3242 const labelList& anchors = cellAnchorPoints_[celli];
3244 writeOBJ(dir, celli, loopPoints(celli), anchors);
Various functions to operate on Lists.
void setSize(const label n)
Alias for resize()
void clear()
Clear the list, i.e. set size to zero.
A HashTable to objects of type <T> with a label key.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
static word timeName(const scalar t, const int precision=precision_)
bool found(const T &val, label pos=0) const
True if the value if found in the list.
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.
label rcIndex(const label i) const noexcept
void size(const label n)
Older name for setAddressableSize.
label fcIndex(const label i) const noexcept
static bool & parRun() noexcept
Test if this a parallel run.
const Cmpt & z() const
Access to the vector z component.
const Cmpt & y() const
Access to the vector y component.
const Cmpt & x() const
Access to the vector x component.
Description of cuts across cells.
const labelListList & cellLoops() const
For each cut cell the cut along the circumference.
void writeCellOBJ(const fileName &dir, const label celli) const
debugging:Write edges of cell and loop
void flipLoopOnly(const label celli)
Flip loop for celli. Does not update anchors. Use with care.
labelList nonAnchorPoints(const labelList &cellPoints, const labelList &anchorPoints, const labelList &loop) const
Invert anchor point selection.
void clearOut()
Clear out demand-driven storage.
Abstract base class. Concrete implementations know how to cut a cell (i.e. determine a loop around th...
Combines edge or vertex in single label. Used to specify cuts across cell circumference.
const polyMesh & mesh() const
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
A face is a list of labels corresponding to mesh vertices.
A class for handling file names.
const Time & time() const
Return the top-level database.
edge cmptType
Component type.
A traits class, which is primarily used for primitives.
void flip()
Reverse the orientation.
Mesh consisting of general polyhedral cells.
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.
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
const labelListList & pointEdges() const
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const labelListList & cellEdges() const
const vectorField & faceCentres() const
label nInternalFaces() const noexcept
Number of internal faces.
const labelListList & cellPoints() const
label nCells() const noexcept
Number of mesh cells.
label nFaces() const noexcept
Number of mesh faces.
const labelListList & pointFaces() const
const labelListList & edgeFaces() const
const labelListList & faceEdges() const
const cellList & cells() const
Patchify triangles based on orientation w.r.t other (triangulated or triangulatable) surfaces.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const labelList nEdges(UPstream::listGatherValues< label >(aMesh.nEdges()))
#define WarningInFunction
Report a warning using Foam::Warning.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
List< label > labelList
A List of labels.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
errorManip< error > abort(error &err)
List< bool > boolList
A List of bools.
static constexpr const zero Zero
Global zero (0)
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
List< edge > edgeList
A List of edges.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< face > faceList
A List of faces.
errorManipArg< error, int > exit(error &err, const int errNo=1)
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]
#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.