54 for (
label i = 0; i < nElems; i++)
75 result[labels[
labelI]] =
true;
101 const Map<label>& map
106 if (!map.found(lst[i]))
117 void Foam::cellCuts::writeUncutOBJ
124 OFstream cutsStream(dir /
"cell_" +
name(cellI) +
".obj");
127 <<
" to " << cutsStream.name() <<
nl;
139 OFstream cutStream(dir /
"cellCuts_" +
name(cellI) +
".obj");
142 <<
" to " << cutStream.name() <<
nl;
148 label pointI = cPoints[i];
149 if (pointIsCut_[pointI])
161 label edgeI = cEdges[i];
163 if (edgeIsCut_[edgeI])
167 const scalar w = edgeWeight_[edgeI];
175 void Foam::cellCuts::writeOBJ
184 OFstream cutsStream(dir /
"cell_" +
name(cellI) +
".obj");
187 <<
" to " << cutsStream.name() <<
nl;
200 OFstream loopStream(dir /
"cellLoop_" +
name(cellI) +
".obj");
203 <<
" to " << loopStream.name() <<
nl;
207 writeOBJ(loopStream, loopPoints, vertI);
211 OFstream anchorStream(dir /
"anchors_" +
name(cellI) +
".obj");
214 <<
" to " << anchorStream.name() <<
endl;
234 label faceI = cFaces[cFaceI];
253 "Foam::cellCuts::edgeEdgeToFace" 254 "(const label cellI, const label edgeA," 255 "const label edgeB) const" 256 ) <<
"cellCuts : Cannot find face on cell " 257 << cellI <<
" that has both edges " << edgeA <<
' ' << edgeB << endl
258 <<
"faces : " << cFaces << endl
259 <<
"edgeA : " <<
mesh().
edges()[edgeA] << endl
260 <<
"edgeB : " <<
mesh().
edges()[edgeB] << endl
261 <<
"Marking the loop across this cell as invalid" <<
endl;
278 label faceI = cFaces[cFaceI];
296 "Foam::cellCuts::edgeVertexToFace" 297 "(const label cellI, const label edgeI, " 298 "const label vertI) const" 299 ) <<
"cellCuts : Cannot find face on cell " 300 << cellI <<
" that has both edge " << edgeI <<
" and vertex " 302 <<
"faces : " << cFaces << endl
303 <<
"edge : " <<
mesh().
edges()[edgeI] << endl
304 <<
"Marking the loop across this cell as invalid" <<
endl;
321 label faceI = cFaces[cFaceI];
335 WarningIn(
"Foam::cellCuts::vertexVertexToFace")
336 <<
"cellCuts : Cannot find face on cell " 337 << cellI <<
" that has vertex " << vertA <<
" and vertex " 339 <<
"faces : " << cFaces << endl
340 <<
"Marking the loop across this cell as invalid" <<
endl;
346 void Foam::cellCuts::calcFaceCuts()
const 362 const face& f = faces[faceI];
386 label vMin1 = f[f.rcIndex(fp)];
391 && !edgeIsCut_[
findEdge(faceI, v0, vMin1)]
394 cuts[cutI++] = vertToEVert(v0);
395 startFp = f.fcIndex(fp);
406 label fp1 = f.fcIndex(fp);
413 if (edgeIsCut_[edgeI] && !pointIsCut_[v0])
415 cuts[cutI++] = edgeToEVert(edgeI);
432 bool allVerticesCut =
true;
436 label fp1 = f.fcIndex(fp);
446 cuts[cutI++] = vertToEVert(v0);
451 allVerticesCut =
false;
454 if (edgeIsCut_[edgeI])
456 cuts[cutI++] = edgeToEVert(edgeI);
465 WarningIn(
"Foam::cellCuts::calcFaceCuts() const")
466 <<
"Face " << faceI <<
" vertices " << f
467 <<
" has all its vertices cut. Not cutting face." <<
endl;
473 if (cutI > 1 && cuts[cutI-1] == cuts[0])
495 const edge& e = edges[fEdges[i]];
499 (e[0] == v0 && e[1] == v1)
500 || (e[0] == v1 && e[1] == v0)
517 const cell& cFaces =
mesh().
cells()[cellI];
521 label faceI = cFaces[cFaceI];
526 bool allOnFace =
true;
534 if (
findIndex(fEdges, getEdge(cut)) == -1)
562 bool Foam::cellCuts::walkPoint
565 const label startCut,
567 const label exclude0,
568 const label exclude1,
570 const label otherCut,
576 label vertI = getVertex(otherCut);
582 label otherFaceI = pFaces[pFaceI];
586 otherFaceI != exclude0
587 && otherFaceI != exclude1
591 label oldNVisited = nVisited;
610 nVisited = oldNVisited;
617 bool Foam::cellCuts::crossEdge
620 const label startCut,
622 const label otherCut,
629 label edgeI = getEdge(otherCut);
634 label oldNVisited = nVisited;
655 nVisited = oldNVisited;
662 bool Foam::cellCuts::addCut
671 if (findPartIndex(visited, nVisited, cut) != -1)
675 truncVisited.setSize(nVisited);
677 Pout<<
"For cell " << cellI <<
" : trying to add duplicate cut " << cut;
679 writeCuts(
Pout, cuts, loopWeights(cuts));
682 writeCuts(
Pout, truncVisited, loopWeights(truncVisited));
689 visited[nVisited++] = cut;
696 bool Foam::cellCuts::walkFace
699 const label startCut,
704 label& beforeLastCut,
709 const labelList& fCuts = faceCuts()[faceI];
711 if (fCuts.size() < 2)
717 if (fCuts.size() == 2)
721 if (!addCut(cellI, cut, nVisited, visited))
733 if (!addCut(cellI, cut, nVisited, visited))
750 for (
label i = 0; i < fCuts.size()-1; i++)
752 if (!addCut(cellI, fCuts[i], nVisited, visited))
757 beforeLastCut = fCuts[fCuts.size()-2];
758 lastCut = fCuts[fCuts.size()-1];
760 else if (fCuts[fCuts.size()-1] == cut)
762 for (
label i = fCuts.size()-1; i >= 1; --i)
764 if (!addCut(cellI, fCuts[i], nVisited, visited))
769 beforeLastCut = fCuts[1];
775 <<
"In middle of cut. cell:" << cellI <<
" face:" << faceI
776 <<
" cuts:" << fCuts <<
" current cut:" << cut <<
endl;
786 bool Foam::cellCuts::walkCell
789 const label startCut,
799 label beforeLastCut = -1;
804 Pout<<
"For cell:" << cellI <<
" walked across face " << faceI
807 writeCuts(
Pout, cuts, loopWeights(cuts));
831 Pout<<
" to last cut ";
833 writeCuts(
Pout, cuts, loopWeights(cuts));
838 if (lastCut == startCut)
846 truncVisited.setSize(nVisited);
848 Pout<<
"For cell " << cellI <<
" : found closed path:";
849 writeCuts(
Pout, truncVisited, loopWeights(truncVisited));
850 Pout<<
" closed by " << lastCut <<
endl;
878 if (isEdge(beforeLastCut))
935 getVertex(beforeLastCut),
979 void Foam::cellCuts::calcCellLoops(
const labelList& cutCells)
994 forAll(allFaceCuts, faceI)
996 const labelList& fCuts = allFaceCuts[faceI];
1003 if (
mesh().isInternalFace(faceI))
1008 else if (fCuts.size() >= 2)
1013 if (
mesh().isInternalFace(faceI))
1027 label cellI = cutCells[i];
1029 bool validLoop =
false;
1032 if (nCutFaces[cellI] >= 3)
1038 Pout<<
"cell:" << cellI <<
" cut faces:" <<
endl;
1041 label faceI = cFaces[i];
1042 const labelList& fCuts = allFaceCuts[faceI];
1044 Pout<<
" face:" << faceI <<
" cuts:";
1045 writeCuts(
Pout, fCuts, loopWeights(fCuts));
1055 label faceI = cFaces[cFaceI];
1057 const labelList& fCuts = allFaceCuts[faceI];
1063 if (fCuts.size() >= 2)
1070 Pout<<
"cell:" << cellI
1071 <<
" start walk at face:" << faceI
1074 writeCuts(
Pout, cuts, loopWeights(cuts));
1109 for (
label i = 0; i < nVisited; i++)
1111 loop[i] = visited[i];
1118 Pout<<
"calcCellLoops(const labelList&) : did not find valid" 1119 <<
" loop for cell " << cellI <<
endl;
1121 writeUncutOBJ(
".", cellI);
1123 cellLoops_[cellI].setSize(0);
1131 cellLoops_[cellI].setSize(0);
1137 void Foam::cellCuts::walkEdges
1143 Map<label>& edgeStatus,
1144 Map<label>& pointStatus
1147 if (pointStatus.insert(pointI, status))
1155 label edgeI = pEdges[pEdgeI];
1160 && edgeStatus.insert(edgeI, status)
1167 walkEdges(cellI, v2, status, edgeStatus, pointStatus);
1186 label pointI = cellPoints[i];
1191 &&
findIndex(loop, vertToEVert(pointI)) == -1
1194 newElems[newElemI++] = pointI;
1204 bool Foam::cellCuts::loopAnchorConsistent
1221 forAll(anchorPoints, ptI)
1225 avg /= anchorPoints.
size();
1228 if (((avg - ctr) & normal) > 0)
1239 bool Foam::cellCuts::calcAnchors
1266 label cut = loop[i];
1270 edgeStatus.insert(getEdge(cut), 0);
1274 pointStatus.insert(getVertex(cut), 0);
1281 label edgeI = cEdges[i];
1282 const edge& e = edges[edgeI];
1284 if (pointStatus.found(e[0]) && pointStatus.found(e[1]))
1286 edgeStatus.insert(edgeI, 0);
1292 label uncutIndex = firstUnique(cPoints, pointStatus);
1294 if (uncutIndex == -1)
1296 WarningIn(
"Foam::cellCuts::calcAnchors")
1297 <<
"Invalid loop " << loop <<
" for cell " << cellI << endl
1298 <<
"Can not find point on cell which is not cut by loop." 1307 walkEdges(cellI, cPoints[uncutIndex], 1, edgeStatus, pointStatus);
1310 uncutIndex = firstUnique(cPoints, pointStatus);
1312 if (uncutIndex == -1)
1316 WarningIn(
"Foam::cellCuts::calcAnchors")
1317 <<
"Invalid loop " << loop <<
" for cell " << cellI << endl
1318 <<
"All vertices of cell are either in loop or in anchor set" 1328 walkEdges(cellI, cPoints[uncutIndex], 2, edgeStatus, pointStatus);
1338 connectedPoints.append(iter.key());
1340 else if (iter() == 2)
1342 otherPoints.append(iter.key());
1345 connectedPoints.shrink();
1346 otherPoints.shrink();
1349 uncutIndex = firstUnique(cPoints, pointStatus);
1351 if (uncutIndex != -1)
1353 WarningIn(
"Foam::cellCuts::calcAnchors")
1354 <<
"Invalid loop " << loop <<
" for cell " << cellI
1355 <<
" since it splits the cell into more than two cells" <<
endl;
1357 writeOBJ(
".", cellI, loopPts, connectedPoints);
1369 label pointI = iter.key();
1379 connectedFaces.insert(pFaces[pFaceI]);
1383 else if (iter() == 2)
1389 otherFaces.insert(pFaces[pFaceI]);
1395 if (connectedFaces.size() < 3)
1397 WarningIn(
"Foam::cellCuts::calcAnchors")
1398 <<
"Invalid loop " << loop <<
" for cell " << cellI
1399 <<
" since would have too few faces on one side." <<
nl 1400 <<
"All faces:" << cFaces <<
endl;
1402 writeOBJ(
".", cellI, loopPts, connectedPoints);
1407 if (otherFaces.size() < 3)
1409 WarningIn(
"Foam::cellCuts::calcAnchors")
1410 <<
"Invalid loop " << loop <<
" for cell " << cellI
1411 <<
" since would have too few faces on one side." <<
nl 1412 <<
"All faces:" << cFaces <<
endl;
1414 writeOBJ(
".", cellI, loopPts, otherPoints);
1427 label faceI = cFaces[i];
1431 bool hasSet1 =
false;
1432 bool hasSet2 =
false;
1434 label prevStat = pointStatus[f[0]];
1439 label pStat = pointStatus[v0];
1441 if (pStat == prevStat)
1444 else if (pStat == 0)
1448 else if (pStat == 1)
1453 WarningIn(
"Foam::cellCuts::calcAnchors")
1454 <<
"Invalid loop " << loop <<
" for cell " << cellI
1455 <<
" since face " << f <<
" would be split into" 1456 <<
" more than two faces" <<
endl;
1458 writeOBJ(
".", cellI, loopPts, otherPoints);
1465 else if (pStat == 2)
1470 WarningIn(
"Foam::cellCuts::calcAnchors")
1471 <<
"Invalid loop " << loop <<
" for cell " << cellI
1472 <<
" since face " << f <<
" would be split into" 1473 <<
" more than two faces" <<
endl;
1475 writeOBJ(
".", cellI, loopPts, otherPoints);
1491 label v1 = f.nextLabel(fp);
1494 label eStat = edgeStatus[edgeI];
1496 if (eStat == prevStat)
1499 else if (eStat == 0)
1503 else if (eStat == 1)
1508 WarningIn(
"Foam::cellCuts::calcAnchors")
1509 <<
"Invalid loop " << loop <<
" for cell " << cellI
1510 <<
" since face " << f <<
" would be split into" 1511 <<
" more than two faces" <<
endl;
1513 writeOBJ(
".", cellI, loopPts, otherPoints);
1520 else if (eStat == 2)
1525 WarningIn(
"Foam::cellCuts::calcAnchors")
1526 <<
"Invalid loop " << loop <<
" for cell " << cellI
1527 <<
" since face " << f <<
" would be split into" 1528 <<
" more than two faces" <<
endl;
1530 writeOBJ(
".", cellI, loopPts, otherPoints);
1546 bool loopOk = loopAnchorConsistent(cellI, loopPts, connectedPoints);
1551 bool otherLoopOk = loopAnchorConsistent(cellI, loopPts, otherPoints);
1553 if (loopOk == otherLoopOk)
1558 WarningIn(
"Foam::cellCuts::calcAnchors")
1559 <<
" For cell:" << cellI
1560 <<
" achorpoints and nonanchorpoints are geometrically" 1561 <<
" on same side!" << endl
1562 <<
"cellPoints:" << cPoints << endl
1563 <<
"loop:" << loop << endl
1564 <<
"anchors:" << connectedPoints << endl
1565 <<
"otherPoints:" << otherPoints <<
endl;
1567 writeOBJ(
".", cellI, loopPts, connectedPoints);
1574 anchorPoints.
transfer(connectedPoints);
1575 connectedPoints.clear();
1579 anchorPoints.
transfer(otherPoints);
1580 otherPoints.clear();
1596 loopPts[fp] = coord(loop[fp], loopWeights[fp]);
1608 label cut = loop[fp];
1612 label edgeI = getEdge(cut);
1614 weights[fp] = edgeWeight_[edgeI];
1618 weights[fp] = -GREAT;
1625 bool Foam::cellCuts::validEdgeLoop
1633 label cut = loop[fp];
1637 label edgeI = getEdge(cut);
1640 if (edgeIsCut_[edgeI])
1646 mag(loopWeights[fp] - edgeWeight_[edgeI])
1677 label vertI = f[fp];
1683 || (
findIndex(loop, vertToEVert(vertI)) != -1)
1695 label edgeI = fEdges[fEdgeI];
1701 || (
findIndex(loop, edgeToEVert(edgeI)) != -1)
1712 bool Foam::cellCuts::conservativeValidLoop
1719 if (loop.
size() < 2)
1726 if (isEdge(loop[cutI]))
1728 label edgeI = getEdge(loop[cutI]);
1730 if (edgeIsCut_[edgeI])
1739 if (pointIsCut_[e.
start()] || pointIsCut_[e.
end()])
1750 label nCuts = countFaceCuts(eFaces[eFaceI], loop);
1763 label vertI = getVertex(loop[cutI]);
1765 if (!pointIsCut_[vertI])
1774 label edgeI = pEdges[pEdgeI];
1776 if (edgeIsCut_[edgeI])
1787 label nCuts = countFaceCuts(pFaces[pFaceI], loop);
1803 bool Foam::cellCuts::validLoop
1818 if (loop.
size() < 2)
1827 if (!conservativeValidLoop(cellI, loop))
1829 Info <<
"Invalid conservative loop: " << loop <<
endl;
1836 label cut = loop[fp];
1837 label nextCut = loop[(fp+1) % loop.
size()];
1840 label meshFaceI = -1;
1844 label edgeI = getEdge(cut);
1848 if (isEdge(nextCut))
1851 label nextEdgeI = getEdge(nextCut);
1854 meshFaceI = edgeEdgeToFace(cellI, edgeI, nextEdgeI);
1856 if (meshFaceI == -1)
1866 label nextVertI = getVertex(nextCut);
1869 if (e.
start() != nextVertI && e.
end() != nextVertI)
1872 meshFaceI = edgeVertexToFace(cellI, edgeI, nextVertI);
1874 if (meshFaceI == -1)
1885 label vertI = getVertex(cut);
1887 if (isEdge(nextCut))
1890 label nextEdgeI = getEdge(nextCut);
1894 if (nextE.
start() != vertI && nextE.
end() != vertI)
1897 meshFaceI = edgeVertexToFace(cellI, nextEdgeI, vertI);
1899 if (meshFaceI == -1)
1908 label nextVertI = getVertex(nextCut);
1913 meshFaceI = vertexVertexToFace(cellI, vertI, nextVertI);
1915 if (meshFaceI == -1)
1923 if (meshFaceI != -1)
1927 edge cutEdge(cut, nextCut);
1931 if (iter == faceSplitCut_.end())
1934 newFaceSplitCut.
insert(meshFaceI, cutEdge);
1939 if (iter() != cutEdge)
1948 label faceContainingLoop = loopFace(cellI, loop);
1950 if (faceContainingLoop != -1)
1953 <<
"Found loop on cell " << cellI <<
" with all points" 1954 <<
" on face " << faceContainingLoop <<
endl;
1967 loopPoints(loop, loopWeights),
1973 void Foam::cellCuts::setFromCellLoops()
1976 pointIsCut_ =
false;
1980 faceSplitCut_.clear();
1982 forAll(cellLoops_, cellI)
1984 const labelList& loop = cellLoops_[cellI];
2009 <<
"Illegal loop " << loop
2010 <<
" when recreating cut-addressing" 2011 <<
" from existing cellLoops for cell " << cellI
2015 cellLoops_[cellI].setSize(0);
2016 cellAnchorPoints_[cellI].setSize(0);
2021 cellAnchorPoints_[cellI].transfer(anchorPoints);
2026 faceSplitCut_.insert(iter.key(), iter());
2032 label cut = loop[cutI];
2036 edgeIsCut_[getEdge(cut)] =
true;
2040 pointIsCut_[getVertex(cut)] =
true;
2048 forAll(edgeIsCut_, edgeI)
2050 if (!edgeIsCut_[edgeI])
2053 edgeWeight_[edgeI] = -GREAT;
2059 bool Foam::cellCuts::setFromCellLoop
2074 str<<
"# edges of cell " << cellI <<
nl;
2088 loopStr<<
"# looppoints for cell " << cellI <<
nl;
2090 pointField pointsOfLoop = loopPoints(loop, loopWeights);
2101 str <<
' ' << i + 1;
2103 str <<
' ' << 1 <<
nl;
2106 bool okLoop =
false;
2108 if (validEdgeLoop(loop, loopWeights))
2116 validLoop(cellI, loop, loopWeights, faceSplitCuts, anchorPoints);
2121 cellLoops_[cellI] = loop;
2122 cellAnchorPoints_[cellI].
transfer(anchorPoints);
2127 faceSplitCut_.insert(iter.key(), iter());
2134 label cut = loop[cutI];
2138 label edgeI = getEdge(cut);
2140 edgeIsCut_[edgeI] =
true;
2142 edgeWeight_[edgeI] = loopWeights[cutI];
2146 label vertI = getVertex(cut);
2148 pointIsCut_[vertI] =
true;
2158 void Foam::cellCuts::setFromCellLoops
2166 pointIsCut_ =
false;
2170 forAll(cellLabels, cellLabelI)
2172 label cellI = cellLabels[cellLabelI];
2174 const labelList& loop = cellLoops[cellLabelI];
2178 const scalarField& loopWeights = cellLoopWeights[cellLabelI];
2180 if (setFromCellLoop(cellI, loop, loopWeights))
2194 void Foam::cellCuts::setFromCellCutter
2201 pointIsCut_ =
false;
2215 forAll(refCells, refCellI)
2217 const refineCell& refCell = refCells[refCellI];
2243 if (setFromCellLoop(cellI, cellLoop, cellLoopWeights))
2249 cellLoops_[cellI].setSize(0);
2251 WarningIn(
"Foam::cellCuts::setFromCellCutter")
2252 <<
"Found loop on cell " << cellI
2253 <<
" that resulted in an unexpected bad cut." 2254 <<
" Suggestions:" <<
nl 2255 <<
" - Turn on the debug switch for 'cellCuts' to get" 2256 <<
" geometry files that identify this cell." <<
nl 2257 <<
" - Also keep in mind to check the defined" 2258 <<
" reference directions, as these are most likely the" 2259 <<
" origin of the problem." 2265 invalidCutCells.
append(cellI);
2266 invalidCutLoops.
append(cellLoop);
2267 invalidCutLoopWeights.
append(cellLoopWeights);
2274 cellLoops_[cellI].setSize(0);
2278 if (debug && invalidCutCells.
size())
2280 invalidCutCells.
shrink();
2281 invalidCutLoops.
shrink();
2282 invalidCutLoopWeights.
shrink();
2284 fileName cutsFile(
"invalidLoopCells.obj");
2286 Pout<<
"cellCuts : writing inValidLoops cells to " << cutsFile <<
endl;
2299 fileName loopsFile(
"invalidLoops.obj");
2301 Pout<<
"cellCuts : writing inValidLoops loops to " << loopsFile <<
endl;
2307 forAll(invalidCutLoops, i)
2312 loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
2320 void Foam::cellCuts::setFromCellCutter
2328 pointIsCut_ =
false;
2344 label cellI = cellLabels[i];
2365 if (setFromCellLoop(cellI, cellLoop, cellLoopWeights))
2371 cellLoops_[cellI].setSize(0);
2376 invalidCutCells.
append(cellI);
2377 invalidCutLoops.
append(cellLoop);
2378 invalidCutLoopWeights.
append(cellLoopWeights);
2385 cellLoops_[cellI].setSize(0);
2389 if (debug && invalidCutCells.
size())
2391 invalidCutCells.
shrink();
2392 invalidCutLoops.
shrink();
2393 invalidCutLoopWeights.
shrink();
2395 fileName cutsFile(
"invalidLoopCells.obj");
2397 Pout<<
"cellCuts : writing inValidLoops cells to " << cutsFile <<
endl;
2410 fileName loopsFile(
"invalidLoops.obj");
2412 Pout<<
"cellCuts : writing inValidLoops loops to " << loopsFile <<
endl;
2418 forAll(invalidCutLoops, i)
2423 loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
2431 void Foam::cellCuts::orientPlanesAndLoops()
2434 forAll(cellLoops_, cellI)
2436 const labelList& loop = cellLoops_[cellI];
2438 if (loop.
size() && cellAnchorPoints_[cellI].empty())
2446 cellAnchorPoints_[cellI]
2453 Pout<<
"cellAnchorPoints:" <<
endl;
2455 forAll(cellAnchorPoints_, cellI)
2457 if (cellLoops_[cellI].size())
2459 if (cellAnchorPoints_[cellI].empty())
2462 <<
"No anchor points for cut cell " << cellI << endl
2468 Pout<<
" cell:" << cellI <<
" anchored at " 2469 << cellAnchorPoints_[cellI] <<
endl;
2477 forAll(cellLoops_, cellI)
2479 if (cellLoops_[cellI].size())
2487 void Foam::cellCuts::calcLoopsAndAddressing(
const labelList& cutCells)
2490 forAll(edgeIsCut_, edgeI)
2492 if (edgeIsCut_[edgeI])
2494 scalar weight = edgeWeight_[edgeI];
2496 if (weight < 0 || weight > 1)
2500 "cellCuts::calcLoopsAndAddressing(const labelList&)" 2501 ) <<
"Weight out of range [0,1]. Edge " << edgeI
2509 edgeWeight_[edgeI] = -GREAT;
2515 calcCellLoops(cutCells);
2520 forAll(cellLoops_, cellI)
2522 const labelList& loop = cellLoops_[cellI];
2526 Pout<<
"cell:" << cellI <<
" ";
2527 writeCuts(
Pout, loop, loopWeights(loop));
2539 void Foam::cellCuts::check()
const 2542 forAll(edgeIsCut_, edgeI)
2544 if (edgeIsCut_[edgeI])
2555 <<
"edge:" << edgeI <<
" vertices:" 2557 <<
" weight:" << edgeWeight_[edgeI] <<
" should have been" 2558 <<
" snapped to one of its endpoints" 2565 if (edgeWeight_[edgeI] > - 1)
2568 <<
"edge:" << edgeI <<
" vertices:" 2570 <<
" weight:" << edgeWeight_[edgeI] <<
" is not cut but" 2571 <<
" its weight is not marked invalid" 2578 forAll(cellLoops_, cellI)
2580 const labelList& loop = cellLoops_[cellI];
2584 label cut = loop[i];
2588 (isEdge(cut) && !edgeIsCut_[getEdge(cut)])
2589 || (!isEdge(cut) && !pointIsCut_[getVertex(cut)])
2593 writeCuts(
Pout, cuts, loopWeights(cuts));
2596 <<
"cell:" << cellI <<
" loop:" 2598 <<
" cut:" << cut <<
" is not marked as cut" 2605 forAll(cellLoops_, cellI)
2607 const labelList& anchors = cellAnchorPoints_[cellI];
2609 const labelList& loop = cellLoops_[cellI];
2614 <<
"cell:" << cellI <<
" loop:" << loop
2615 <<
" has no anchor points" 2622 label cut = loop[i];
2627 &&
findIndex(anchors, getVertex(cut)) != -1
2631 <<
"cell:" << cellI <<
" loop:" << loop
2632 <<
" anchor points:" << anchors
2633 <<
" anchor:" << getVertex(cut) <<
" is part of loop" 2643 label faceI = iter.key();
2645 if (
mesh().isInternalFace(faceI))
2650 if (cellLoops_[own].empty() && cellLoops_[nei].empty())
2653 <<
"Internal face:" << faceI <<
" cut by " << iter()
2654 <<
" has owner:" << own
2655 <<
" and neighbour:" << nei
2656 <<
" that are both uncut" 2664 if (cellLoops_[own].empty())
2667 <<
"Boundary face:" << faceI <<
" cut by " << iter()
2668 <<
" has owner:" << own
2679 Foam::cellCuts::cellCuts
2691 edgeWeight_(
expand(mesh.
nEdges(), meshEdges, meshEdgeWeights)),
2693 faceSplitCut_(cutCells.
size()),
2694 cellLoops_(mesh.
nCells()),
2696 cellAnchorPoints_(mesh.
nCells())
2700 Pout<<
"cellCuts : constructor from cut verts and edges" <<
endl;
2703 calcLoopsAndAddressing(cutCells);
2706 orientPlanesAndLoops();
2717 Pout<<
"cellCuts : leaving constructor from cut verts and edges" 2723 Foam::cellCuts::cellCuts
2734 edgeWeight_(
expand(mesh.
nEdges(), meshEdges, meshEdgeWeights)),
2736 faceSplitCut_(mesh.
nFaces()/10 + 1),
2737 cellLoops_(mesh.
nCells()),
2739 cellAnchorPoints_(mesh.
nCells())
2746 Pout<<
"cellCuts : constructor from cellLoops" <<
endl;
2752 orientPlanesAndLoops();
2763 Pout<<
"cellCuts : leaving constructor from cellLoops" <<
endl;
2768 Foam::cellCuts::cellCuts
2777 pointIsCut_(mesh.
nPoints(),
false),
2778 edgeIsCut_(mesh.
nEdges(),
false),
2779 edgeWeight_(mesh.
nEdges(), -GREAT),
2781 faceSplitCut_(cellLabels.
size()),
2782 cellLoops_(mesh.
nCells()),
2784 cellAnchorPoints_(mesh.
nCells())
2788 Pout<<
"cellCuts : constructor from cellLoops" <<
endl;
2793 setFromCellLoops(cellLabels, cellLoops, cellEdgeWeights);
2796 orientPlanesAndLoops();
2807 Pout<<
"cellCuts : leaving constructor from cellLoops" <<
endl;
2812 Foam::cellCuts::cellCuts
2820 pointIsCut_(mesh.
nPoints(),
false),
2821 edgeIsCut_(mesh.
nEdges(),
false),
2822 edgeWeight_(mesh.
nEdges(), -GREAT),
2824 faceSplitCut_(refCells.
size()),
2825 cellLoops_(mesh.
nCells()),
2827 cellAnchorPoints_(mesh.
nCells())
2831 Pout<<
"cellCuts : constructor from cellCutter" <<
endl;
2836 setFromCellCutter(cellCutter, refCells);
2839 orientPlanesAndLoops();
2850 Pout<<
"cellCuts : leaving constructor from cellCutter" <<
endl;
2855 Foam::cellCuts::cellCuts
2864 pointIsCut_(mesh.
nPoints(),
false),
2865 edgeIsCut_(mesh.
nEdges(),
false),
2866 edgeWeight_(mesh.
nEdges(), -GREAT),
2868 faceSplitCut_(cellLabels.
size()),
2869 cellLoops_(mesh.
nCells()),
2871 cellAnchorPoints_(mesh.
nCells())
2875 Pout<<
"cellCuts : constructor from cellCutter with prescribed plane" 2881 setFromCellCutter(cellCutter, cellLabels, cutPlanes);
2884 orientPlanesAndLoops();
2895 Pout<<
"cellCuts : leaving constructor from cellCutter with prescribed" 2896 <<
" plane" <<
endl;
2901 Foam::cellCuts::cellCuts
2914 pointIsCut_(pointIsCut),
2915 edgeIsCut_(edgeIsCut),
2916 edgeWeight_(edgeWeight),
2918 faceSplitCut_(faceSplitCut),
2919 cellLoops_(cellLoops),
2921 cellAnchorPoints_(cellAnchorPoints)
2925 Pout<<
"cellCuts : constructor from components" <<
endl;
2926 Pout<<
"cellCuts : leaving constructor from components" <<
endl;
2949 const labelList& loop = cellLoops_[cellI];
2955 label cut = loop[fp];
2959 label edgeI = getEdge(cut);
2961 loopPts[fp] = coord(cut, edgeWeight_[edgeI]);
2965 loopPts[fp] = coord(cut, -GREAT);
2979 cellAnchorPoints_[cellI] =
2982 mesh().cellPoints()[cellI],
2983 cellAnchorPoints_[cellI],
2997 void Foam::cellCuts::writeOBJ
3004 label startVertI = vertI;
3008 const point& pt = loopPts[fp];
3010 os <<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z() <<
endl;
3018 os <<
' ' << startVertI + fp + 1;
3028 forAll(cellLoops_, cellI)
3030 writeOBJ(os, loopPoints(cellI), vertI);
3037 const labelList& anchors = cellAnchorPoints_[cellI];
3039 writeOBJ(dir, cellI, loopPoints(cellI), anchors);
const labelListList & cellPoints() const
virtual bool cut(const vector &refDir, const label cellI, const boolList &vertIsCut, const boolList &edgeIsCut, const scalarField &edgeWeight, labelList &loop, scalarField &loopWeights) const =0
Create cut along circumference of cellI. Gets current mesh cuts.
const pointField & points
Container with cells to refine. Refinement given as single direction.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
void clearOut()
Clear out demand driven storage.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
dimensioned< scalar > mag(const dimensioned< Type > &)
An STL-conforming const_iterator.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
word name(const complex &)
Return a string representation of a complex.
A cell is defined as a list of faces with extra functionality.
const labelListList & pointEdges() const
bool empty() const
Return true if the UList is empty (ie, size() is zero).
bool insert(const label &, const T &newElmt)
Insert a new hashedEntry.
void getEdgeFaces(const primitiveMesh &, const label cellI, const label edgeI, label &face0, label &face1)
Get faces on cell using edgeI. Throws error if no two found.
void deleteDemandDrivenData(DataPtr &dataPtr)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
void size(const label)
Override size to be inconsistent with allocated storage.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
const cellList & cells() const
Various functions to operate on Lists.
void writeCellOBJ(const fileName &dir, const label cellI) const
debugging:Write edges of cell and loop
bool faceOnCell(const primitiveMesh &, const label cellI, const label faceI)
Is face used by cell.
static const labelSphericalTensor labelI(1)
Identity labelTensor.
void flipLoopOnly(const label cellI)
Flip loop for cellI. Does not update anchors. Use with care.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
point centre(const pointField &) const
Centre point of face.
const labelListList & faceEdges() const
vectorField pointField
pointField is a vectorField.
const labelListList & cellEdges() const
string expand(const string &, const HashTable< string, word, string::hash > &mapping, const char sigil= '$')
Expand occurences of variables according to the mapping.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const Time & time() const
Return the top-level database.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
bool edgeOnCell(const primitiveMesh &, const label cellI, const label edgeI)
Is edge used by cell.
A face is a list of labels corresponding to mesh vertices.
void flip(const label cellI)
Flip loop for cellI. Updates anchor points as well.
void setSize(const label)
Reset size of List.
Ostream & endl(Ostream &os)
Add newline and flush stream.
#define WarningIn(functionName)
Report a warning using Foam::Warning.
virtual const pointField & points() const
Return raw points.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
static const label labelMin
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
const labelListList & pointFaces() const
errorManip< error > abort(error &err)
virtual const labelList & faceOwner() const
Return face owner.
label otherFace(const primitiveMesh &, const label cellI, const label faceI, const label edgeI)
Return face on cell using edgeI but not faceI. Throws error.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
const vector & direction() const
const cellShapeList & cells
Abstract base class. Concrete implementations know how to cut a cell (i.e. determine a loop around th...
Mesh consisting of general polyhedral cells.
vector normal(const pointField &) const
Vector normal; magnitude is equal to area of face.
const labelListList & edgeFaces() const
label end() const
Return end vertex label.
List< label > labelList
A List of labels.
A class for handling file names.
A normal distribution model.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
void reverse(UList< T > &, const label n)
label walkFace(const primitiveMesh &, const label faceI, const label startEdgeI, const label startVertI, const label nEdges)
Returns label of edge nEdges away from startEdge (in the direction.
label start() const
Return start vertex label.
virtual const faceList & faces() const
Return raw faces.
label findEdge(const edgeList &edges, const labelList &candidates, const label v0, const label v1)
Return edge among candidates that uses the two vertices.
List< labelList > labelListList
A List of labelList.
virtual const labelList & faceNeighbour() const
Return face neighbour.
labelList nonAnchorPoints(const labelList &cellPoints, const labelList &anchorPoints, const labelList &loop) const
Invert anchor point selection.
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Combines edge or vertex in single label. Used to specify cuts across cell circumference.
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
List< bool > boolList
Bool container classes.