63 Foam::label Foam::cellCuts::findPartIndex
70 for (label i = 0; i < nElems; i++)
91 result[labels[labelI]] =
true;
108 result[labels[labelI]] = weights[labelI];
114 Foam::label Foam::cellCuts::firstUnique
117 const Map<label>& map
122 if (!map.found(lst[i]))
133 void 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 " 280 void 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];
338 void 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;
386 Foam::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;
422 Foam::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;
457 Foam::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;
489 void 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])
627 Foam::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)
656 Foam::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)))
688 if (!
f.
found(getVertex(cut)))
707 bool 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;
762 bool Foam::cellCuts::crossEdge
765 const label startCut,
767 const label otherCut,
774 label edgeI =
getEdge(otherCut);
779 label oldNVisited = nVisited;
800 nVisited = oldNVisited;
807 bool 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;
845 bool 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;
938 bool 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),
1131 void Foam::cellCuts::calcCellLoops(
const labelList& cutCells)
1146 forAll(allFaceCuts, facei)
1148 const labelList& fCuts = allFaceCuts[facei];
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();
1292 void 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;
1359 bool 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)
1392 bool 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];
1795 weights[fp] = edgeWeight_[edgeI];
1799 weights[fp] = -GREAT;
1806 bool Foam::cellCuts::validEdgeLoop
1814 label cut = loop[fp];
1821 if (edgeIsCut_[edgeI])
1827 mag(loopWeights[fp] - edgeWeight_[edgeI])
1841 Foam::label Foam::cellCuts::countFaceCuts
1858 label vertI =
f[fp];
1864 || loop.found(vertToEVert(vertI))
1876 label edgeI = fEdges[fEdgeI];
1882 || loop.found(edgeToEVert(edgeI))
1893 bool 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);
1984 bool 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;
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),
2157 void 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;
2243 bool 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;
2342 void 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();
2378 void 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]),
2504 void 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]),
2615 void 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())
2671 void 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));
2721 void 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
2871 Foam::cellCuts::cellCuts
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" 2916 Foam::cellCuts::cellCuts
2928 edgeIsCut_(
expand(
mesh.nEdges(), meshEdges)),
2929 edgeWeight_(
expand(
mesh.nEdges(), meshEdges, meshEdgeWeights)),
2930 faceSplitCut_(
mesh.nFaces()/10 + 1),
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;
2965 Foam::cellCuts::cellCuts
2977 edgeIsCut_(
mesh.nEdges(), false),
2978 edgeWeight_(
mesh.nEdges(), -GREAT),
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;
3013 Foam::cellCuts::cellCuts
3024 edgeIsCut_(
mesh.nEdges(), false),
3025 edgeWeight_(
mesh.nEdges(), -GREAT),
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;
3060 Foam::cellCuts::cellCuts
3072 edgeIsCut_(
mesh.nEdges(), false),
3073 edgeWeight_(
mesh.nEdges(), -GREAT),
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;
3110 Foam::cellCuts::cellCuts
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];
3163 loopPts[fp] = coord(cut, edgeWeight_[edgeI]);
3167 loopPts[fp] = coord(cut, -GREAT);
3181 cellAnchorPoints_[celli] =
3184 mesh().cellPoints()[celli],
3185 cellAnchorPoints_[celli],
3199 void 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;
3226 void Foam::cellCuts::writeOBJ(
Ostream&
os)
const 3230 forAll(cellLoops_, celli)
3232 if (cellLoops_[celli].size())
3242 const labelList& anchors = cellAnchorPoints_[celli];
3244 writeOBJ(dir, celli, loopPoints(celli), anchors);
const polyBoundaryMesh & pbm
const labelListList & cellEdges() const
void size(const label n)
Older name for setAddressableSize.
A class for handling file names.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
const labelListList & faceEdges() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual const labelList & faceNeighbour() const
Return face neighbour.
List< edge > edgeList
List of edge.
const labelListList & pointEdges() const
constexpr char nl
The newline '\n' character (0x0a)
void flip(const label celli)
Flip loop for celli. Updates anchor points as well.
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool found(const T &val, label pos=0) const
Same as contains()
Abstract base class. Concrete implementations know how to cut a cell (i.e. determine a loop around th...
A traits class, which is primarily used for primitives and vector-space.
static bool & parRun() noexcept
Test if this a parallel run.
const cellList & cells() const
const Time & time() const
Return the top-level database.
label nFaces() const noexcept
Number of mesh faces.
List< labelList > labelListList
List of labelList.
Various functions to operate on Lists.
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...
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
List< face > faceList
List of faces.
static label getEdge(List< DynamicList< label >> &pe, DynamicList< edge > &es, const label pointi, const label nextPointi)
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
An edge is a list of two vertex labels. This can correspond to a directed graph edge or an edge on a ...
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
void clear()
Clear the list, i.e. set size to zero.
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void clearOut()
Clear out demand-driven storage.
virtual const labelList & faceOwner() const
Return face owner.
Combines edge or vertex in single label. Used to specify cuts across cell circumference.
label nInternalFaces() const noexcept
Number of internal faces.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
virtual const faceList & faces() const
Return raw faces.
errorManip< error > abort(error &err)
label find(const T &val) const
Find index of the first occurrence of the value.
const polyMesh & mesh() const
An Ostream is an abstract base class for all output systems (streams, files, token lists...
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
static word timeName(const scalar t, const int precision=precision_)
Return a time name for the given scalar time value formatted with the given precision.
int debug
Static debugging option.
OBJstream os(runTime.globalPath()/outputName)
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]
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const vectorField & faceCentres() const
const labelListList & cellLoops() const
For each cut cell the cut along the circumference.
#define WarningInFunction
Report a warning using Foam::Warning.
label nCells() const noexcept
Number of mesh cells.
messageStream Info
Information stream (stdout output on master, null elsewhere)
const labelListList & pointFaces() const
Mesh consisting of general polyhedral cells.
List< label > labelList
A List of labels.
label rcIndex(const label i) const noexcept
The reverse circular index. The previous index in the list which returns to the last at the beginning...
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
List< bool > boolList
A List of bools.
void writeCellOBJ(const fileName &dir, const label celli) const
debugging:Write edges of cell and loop
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
const labelListList & edgeFaces() const
const labelListList & cellPoints() const
labelList nonAnchorPoints(const labelList &cellPoints, const labelList &anchorPoints, const labelList &loop) const
Invert anchor point selection.
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
forAllConstIters(mixture.phases(), phase)
void flipLoopOnly(const label celli)
Flip loop for celli. Does not update anchors. Use with care.
A HashTable to objects of type <T> with a label key.
string expand(const std::string &s, const HashTable< string > &mapping, const char sigil='$')
Expand occurrences of variables according to the mapping and return the expanded string.
static constexpr const zero Zero
Global zero (0)