45 Foam::label Foam::cellCuts::findPartIndex
52 for (label i = 0; i < nElems; i++)
73 result[labels[
labelI]] =
true;
97 Foam::label Foam::cellCuts::firstUnique
100 const Map<label>& map
105 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];
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;
224 Foam::label Foam::cellCuts::edgeEdgeToFace
235 label faceI = cFaces[cFaceI];
254 "Foam::cellCuts::edgeEdgeToFace"
255 "(const label cellI, const label edgeA,"
256 "const label edgeB) const"
257 ) <<
"cellCuts : Cannot find face on cell "
258 << cellI <<
" that has both edges " << edgeA <<
' ' << edgeB <<
endl
259 <<
"faces : " << cFaces <<
endl
262 <<
"Marking the loop across this cell as invalid" <<
endl;
269 Foam::label Foam::cellCuts::edgeVertexToFace
280 label faceI = cFaces[cFaceI];
298 "Foam::cellCuts::edgeVertexToFace"
299 "(const label cellI, const label edgeI, "
300 "const label vertI) const"
301 ) <<
"cellCuts : Cannot find face on cell "
302 << cellI <<
" that has both edge " << edgeI <<
" and vertex "
304 <<
"faces : " << cFaces <<
endl
306 <<
"Marking the loop across this cell as invalid" <<
endl;
313 Foam::label Foam::cellCuts::vertexVertexToFace
324 label faceI = cFaces[cFaceI];
338 WarningIn(
"Foam::cellCuts::vertexVertexToFace")
339 <<
"cellCuts : Cannot find face on cell "
340 << cellI <<
" that has vertex " << vertA <<
" and vertex "
342 <<
"faces : " << cFaces <<
endl
343 <<
"Marking the loop across this cell as invalid" <<
endl;
349 void Foam::cellCuts::calcFaceCuts()
const
363 for (label faceI = 0; faceI <
mesh().
nFaces(); faceI++)
365 const face& f = faces[faceI];
389 label vMin1 = f[f.rcIndex(fp)];
394 && !edgeIsCut_[
findEdge(faceI, v0, vMin1)]
397 cuts[cutI++] = vertToEVert(v0);
398 startFp = f.fcIndex(fp);
409 label fp1 = f.fcIndex(fp);
414 label edgeI =
findEdge(faceI, v0, v1);
416 if (edgeIsCut_[edgeI] && !pointIsCut_[v0])
418 cuts[cutI++] = edgeToEVert(edgeI);
435 bool allVerticesCut =
true;
439 label fp1 = f.fcIndex(fp);
445 label edgeI =
findEdge(faceI, v0, v1);
449 cuts[cutI++] = vertToEVert(v0);
454 allVerticesCut =
false;
457 if (edgeIsCut_[edgeI])
459 cuts[cutI++] = edgeToEVert(edgeI);
468 WarningIn(
"Foam::cellCuts::calcFaceCuts() const")
469 <<
"Face " << faceI <<
" vertices " << f
470 <<
" has all its vertices cut. Not cutting face." <<
endl;
476 if (cutI > 1 && cuts[cutI-1] == cuts[0])
499 const edge&
e = edges[fEdges[i]];
503 (e[0] == v0 && e[1] == v1)
504 || (e[0] == v1 && e[1] == v0)
516 Foam::label Foam::cellCuts::loopFace
522 const cell& cFaces =
mesh().
cells()[cellI];
526 label faceI = cFaces[cFaceI];
531 bool allOnFace =
true;
539 if (
findIndex(fEdges, getEdge(cut)) == -1)
568 bool Foam::cellCuts::walkPoint
571 const label startCut,
573 const label exclude0,
574 const label exclude1,
576 const label otherCut,
582 label vertI = getVertex(otherCut);
588 label otherFaceI = pFaces[pFaceI];
592 otherFaceI != exclude0
593 && otherFaceI != exclude1
597 label oldNVisited = nVisited;
616 nVisited = oldNVisited;
624 bool Foam::cellCuts::crossEdge
627 const label startCut,
629 const label otherCut,
636 label edgeI = getEdge(otherCut);
641 label oldNVisited = nVisited;
662 nVisited = oldNVisited;
669 bool Foam::cellCuts::addCut
678 if (findPartIndex(visited, nVisited, cut) != -1)
682 truncVisited.setSize(nVisited);
684 Pout<<
"For cell " << cellI <<
" : trying to add duplicate cut " << cut;
686 writeCuts(
Pout, cuts, loopWeights(cuts));
689 writeCuts(
Pout, truncVisited, loopWeights(truncVisited));
696 visited[nVisited++] = cut;
708 const label startCut,
713 label& beforeLastCut,
718 const labelList& fCuts = faceCuts()[faceI];
720 if (fCuts.size() < 2)
726 if (fCuts.size() == 2)
730 if (!addCut(cellI, cut, nVisited, visited))
742 if (!addCut(cellI, cut, nVisited, visited))
759 for (label i = 0; i < fCuts.size()-1; i++)
761 if (!addCut(cellI, fCuts[i], nVisited, visited))
766 beforeLastCut = fCuts[fCuts.size()-2];
767 lastCut = fCuts[fCuts.size()-1];
769 else if (fCuts[fCuts.size()-1] == cut)
771 for (label i = fCuts.size()-1; i >= 1; --i)
773 if (!addCut(cellI, fCuts[i], nVisited, visited))
778 beforeLastCut = fCuts[1];
784 <<
"In middle of cut. cell:" << cellI <<
" face:" << faceI
785 <<
" cuts:" << fCuts <<
" current cut:" << cut <<
endl;
797 bool Foam::cellCuts::walkCell
800 const label startCut,
810 label beforeLastCut = -1;
815 Pout<<
"For cell:" << cellI <<
" walked across face " << faceI
818 writeCuts(
Pout, cuts, loopWeights(cuts));
842 Pout<<
" to last cut ";
844 writeCuts(
Pout, cuts, loopWeights(cuts));
849 if (lastCut == startCut)
857 truncVisited.setSize(nVisited);
859 Pout<<
"For cell " << cellI <<
" : found closed path:";
860 writeCuts(
Pout, truncVisited, loopWeights(truncVisited));
861 Pout<<
" closed by " << lastCut <<
endl;
889 if (isEdge(beforeLastCut))
946 getVertex(beforeLastCut),
995 void Foam::cellCuts::calcCellLoops(
const labelList& cutCells)
1004 forAll(allFaceCuts, faceI)
1006 const labelList& fCuts = allFaceCuts[faceI];
1013 if (
mesh().isInternalFace(faceI))
1018 else if (fCuts.size() >= 2)
1023 if (
mesh().isInternalFace(faceI))
1037 label cellI = cutCells[i];
1039 bool validLoop =
false;
1042 if (nCutFaces[cellI] >= 3)
1048 Pout<<
"cell:" << cellI <<
" cut faces:" <<
endl;
1051 label faceI = cFaces[i];
1052 const labelList& fCuts = allFaceCuts[faceI];
1054 Pout<<
" face:" << faceI <<
" cuts:";
1055 writeCuts(
Pout, fCuts, loopWeights(fCuts));
1065 label faceI = cFaces[cFaceI];
1067 const labelList& fCuts = allFaceCuts[faceI];
1073 if (fCuts.size() >= 2)
1080 Pout<<
"cell:" << cellI
1081 <<
" start walk at face:" << faceI
1084 writeCuts(
Pout, cuts, loopWeights(cuts));
1119 for (label i = 0; i < nVisited; i++)
1121 loop[i] = visited[i];
1128 Pout<<
"calcCellLoops(const labelList&) : did not find valid"
1129 <<
" loop for cell " << cellI <<
endl;
1131 writeUncutOBJ(
".", cellI);
1133 cellLoops_[cellI].setSize(0);
1141 cellLoops_[cellI].setSize(0);
1149 void Foam::cellCuts::walkEdges
1155 Map<label>& edgeStatus,
1156 Map<label>& pointStatus
1159 if (pointStatus.insert(pointI, status))
1167 label edgeI = pEdges[pEdgeI];
1172 && edgeStatus.insert(edgeI, status)
1177 label
v2 =
mesh().
edges()[edgeI].otherVertex(pointI);
1179 walkEdges(cellI, v2, status, edgeStatus, pointStatus);
1199 label pointI = cellPoints[i];
1204 &&
findIndex(loop, vertToEVert(pointI)) == -1
1207 newElems[newElemI++] = pointI;
1218 bool Foam::cellCuts::loopAnchorConsistent
1235 forAll(anchorPoints, ptI)
1239 avg /= anchorPoints.
size();
1242 if (((avg - ctr) & normal) > 0)
1257 bool Foam::cellCuts::calcAnchors
1270 const cell& cFaces =
mesh().
cells()[cellI];
1278 Map<label> pointStatus(2*cPoints.size());
1279 Map<label> edgeStatus(2*cEdges.size());
1284 label cut = loop[i];
1288 edgeStatus.insert(getEdge(cut), 0);
1292 pointStatus.insert(getVertex(cut), 0);
1299 label edgeI = cEdges[i];
1300 const edge& e = edges[edgeI];
1302 if (pointStatus.found(e[0]) && pointStatus.found(e[1]))
1304 edgeStatus.insert(edgeI, 0);
1310 label uncutIndex = firstUnique(cPoints, pointStatus);
1312 if (uncutIndex == -1)
1314 WarningIn(
"Foam::cellCuts::calcAnchors")
1315 <<
"Invalid loop " << loop <<
" for cell " << cellI << endl
1316 <<
"Can not find point on cell which is not cut by loop."
1325 walkEdges(cellI, cPoints[uncutIndex], 1, edgeStatus, pointStatus);
1328 uncutIndex = firstUnique(cPoints, pointStatus);
1330 if (uncutIndex == -1)
1334 WarningIn(
"Foam::cellCuts::calcAnchors")
1335 <<
"Invalid loop " << loop <<
" for cell " << cellI << endl
1336 <<
"All vertices of cell are either in loop or in anchor set"
1346 walkEdges(cellI, cPoints[uncutIndex], 2, edgeStatus, pointStatus);
1349 DynamicList<label> connectedPoints(cPoints.size());
1350 DynamicList<label> otherPoints(cPoints.size());
1356 connectedPoints.append(iter.key());
1358 else if (iter() == 2)
1360 otherPoints.append(iter.key());
1363 connectedPoints.shrink();
1364 otherPoints.shrink();
1367 uncutIndex = firstUnique(cPoints, pointStatus);
1369 if (uncutIndex != -1)
1371 WarningIn(
"Foam::cellCuts::calcAnchors")
1372 <<
"Invalid loop " << loop <<
" for cell " << cellI
1373 <<
" since it splits the cell into more than two cells" <<
endl;
1375 writeOBJ(
".", cellI, loopPts, connectedPoints);
1387 label pointI = iter.key();
1397 connectedFaces.insert(pFaces[pFaceI]);
1401 else if (iter() == 2)
1407 otherFaces.insert(pFaces[pFaceI]);
1413 if (connectedFaces.size() < 3)
1415 WarningIn(
"Foam::cellCuts::calcAnchors")
1416 <<
"Invalid loop " << loop <<
" for cell " << cellI
1417 <<
" since would have too few faces on one side." <<
nl
1418 <<
"All faces:" << cFaces <<
endl;
1420 writeOBJ(
".", cellI, loopPts, connectedPoints);
1425 if (otherFaces.size() < 3)
1427 WarningIn(
"Foam::cellCuts::calcAnchors")
1428 <<
"Invalid loop " << loop <<
" for cell " << cellI
1429 <<
" since would have too few faces on one side." <<
nl
1430 <<
"All faces:" << cFaces <<
endl;
1432 writeOBJ(
".", cellI, loopPts, otherPoints);
1445 label faceI = cFaces[i];
1449 bool hasSet1 =
false;
1450 bool hasSet2 =
false;
1452 label prevStat = pointStatus[f[0]];
1457 label pStat = pointStatus[v0];
1459 if (pStat == prevStat)
1462 else if (pStat == 0)
1466 else if (pStat == 1)
1471 WarningIn(
"Foam::cellCuts::calcAnchors")
1472 <<
"Invalid loop " << loop <<
" for cell " << cellI
1473 <<
" since face " << f <<
" would be split into"
1474 <<
" more than two faces" <<
endl;
1476 writeOBJ(
".", cellI, loopPts, otherPoints);
1483 else if (pStat == 2)
1488 WarningIn(
"Foam::cellCuts::calcAnchors")
1489 <<
"Invalid loop " << loop <<
" for cell " << cellI
1490 <<
" since face " << f <<
" would be split into"
1491 <<
" more than two faces" <<
endl;
1493 writeOBJ(
".", cellI, loopPts, otherPoints);
1509 label v1 = f.nextLabel(fp);
1510 label edgeI =
findEdge(faceI, v0, v1);
1512 label eStat = edgeStatus[edgeI];
1514 if (eStat == prevStat)
1517 else if (eStat == 0)
1521 else if (eStat == 1)
1526 WarningIn(
"Foam::cellCuts::calcAnchors")
1527 <<
"Invalid loop " << loop <<
" for cell " << cellI
1528 <<
" since face " << f <<
" would be split into"
1529 <<
" more than two faces" <<
endl;
1531 writeOBJ(
".", cellI, loopPts, otherPoints);
1538 else if (eStat == 2)
1543 WarningIn(
"Foam::cellCuts::calcAnchors")
1544 <<
"Invalid loop " << loop <<
" for cell " << cellI
1545 <<
" since face " << f <<
" would be split into"
1546 <<
" more than two faces" <<
endl;
1548 writeOBJ(
".", cellI, loopPts, otherPoints);
1564 bool loopOk = loopAnchorConsistent(cellI, loopPts, connectedPoints);
1569 bool otherLoopOk = loopAnchorConsistent(cellI, loopPts, otherPoints);
1571 if (loopOk == otherLoopOk)
1576 WarningIn(
"Foam::cellCuts::calcAnchors")
1577 <<
" For cell:" << cellI
1578 <<
" achorpoints and nonanchorpoints are geometrically"
1579 <<
" on same side!" << endl
1580 <<
"cellPoints:" << cPoints << endl
1581 <<
"loop:" << loop << endl
1582 <<
"anchors:" << connectedPoints << endl
1583 <<
"otherPoints:" << otherPoints <<
endl;
1585 writeOBJ(
".", cellI, loopPts, connectedPoints);
1592 anchorPoints.transfer(connectedPoints);
1593 connectedPoints.clear();
1597 anchorPoints.transfer(otherPoints);
1598 otherPoints.clear();
1614 loopPts[fp] = coord(loop[fp], loopWeights[fp]);
1627 label cut = loop[fp];
1631 label edgeI = getEdge(cut);
1633 weights[fp] = edgeWeight_[edgeI];
1637 weights[fp] = -GREAT;
1645 bool Foam::cellCuts::validEdgeLoop
1653 label cut = loop[fp];
1657 label edgeI = getEdge(cut);
1660 if (edgeIsCut_[edgeI])
1666 mag(loopWeights[fp] - edgeWeight_[edgeI])
1683 Foam::label Foam::cellCuts::countFaceCuts
1696 label vertI = f[fp];
1702 || (
findIndex(loop, vertToEVert(vertI)) != -1)
1714 label edgeI = fEdges[fEdgeI];
1720 || (
findIndex(loop, edgeToEVert(edgeI)) != -1)
1733 bool Foam::cellCuts::conservativeValidLoop
1740 if (loop.size() < 2)
1747 if (isEdge(loop[cutI]))
1749 label edgeI = getEdge(loop[cutI]);
1751 if (edgeIsCut_[edgeI])
1760 if (pointIsCut_[e.start()] || pointIsCut_[e.end()])
1771 label nCuts = countFaceCuts(eFaces[eFaceI], loop);
1784 label vertI = getVertex(loop[cutI]);
1786 if (!pointIsCut_[vertI])
1795 label edgeI = pEdges[pEdgeI];
1797 if (edgeIsCut_[edgeI])
1808 label nCuts = countFaceCuts(pFaces[pFaceI], loop);
1828 bool Foam::cellCuts::validLoop
1834 Map<edge>& newFaceSplitCut,
1838 if (loop.size() < 2)
1847 if (!conservativeValidLoop(cellI, loop))
1855 label cut = loop[fp];
1856 label nextCut = loop[(fp+1) % loop.size()];
1859 label meshFaceI = -1;
1863 label edgeI = getEdge(cut);
1867 if (isEdge(nextCut))
1870 label nextEdgeI = getEdge(nextCut);
1873 meshFaceI = edgeEdgeToFace(cellI, edgeI, nextEdgeI);
1875 if (meshFaceI == -1)
1885 label nextVertI = getVertex(nextCut);
1888 if (e.start() != nextVertI && e.end() != nextVertI)
1891 meshFaceI = edgeVertexToFace(cellI, edgeI, nextVertI);
1893 if (meshFaceI == -1)
1904 label vertI = getVertex(cut);
1906 if (isEdge(nextCut))
1909 label nextEdgeI = getEdge(nextCut);
1911 const edge& nextE =
mesh().
edges()[nextEdgeI];
1913 if (nextE.start() != vertI && nextE.end() != vertI)
1916 meshFaceI = edgeVertexToFace(cellI, nextEdgeI, vertI);
1918 if (meshFaceI == -1)
1927 label nextVertI = getVertex(nextCut);
1932 meshFaceI = vertexVertexToFace(cellI, vertI, nextVertI);
1934 if (meshFaceI == -1)
1942 if (meshFaceI != -1)
1946 edge cutEdge(cut, nextCut);
1950 if (iter == faceSplitCut_.end())
1953 newFaceSplitCut.insert(meshFaceI, cutEdge);
1958 if (iter() != cutEdge)
1967 label faceContainingLoop = loopFace(cellI, loop);
1969 if (faceContainingLoop != -1)
1972 <<
"Found loop on cell " << cellI <<
" with all points"
1973 <<
" on face " << faceContainingLoop <<
endl;
1986 loopPoints(loop, loopWeights),
1995 void Foam::cellCuts::setFromCellLoops()
1998 pointIsCut_ =
false;
2002 faceSplitCut_.clear();
2004 forAll(cellLoops_, cellI)
2006 const labelList& loop = cellLoops_[cellI];
2011 Map<edge> faceSplitCuts(loop.size());
2031 <<
"Illegal loop " << loop
2032 <<
" when recreating cut-addressing"
2033 <<
" from existing cellLoops for cell " << cellI
2037 cellLoops_[cellI].setSize(0);
2038 cellAnchorPoints_[cellI].setSize(0);
2043 cellAnchorPoints_[cellI].transfer(anchorPoints);
2048 faceSplitCut_.insert(iter.key(), iter());
2054 label cut = loop[cutI];
2058 edgeIsCut_[getEdge(cut)] =
true;
2062 pointIsCut_[getVertex(cut)] =
true;
2070 forAll(edgeIsCut_, edgeI)
2072 if (!edgeIsCut_[edgeI])
2075 edgeWeight_[edgeI] = -GREAT;
2083 bool Foam::cellCuts::setFromCellLoop
2093 OFstream str(
"last_cell.obj");
2095 str<<
"# edges of cell " << cellI <<
nl;
2107 OFstream loopStr(
"last_loop.obj");
2109 loopStr<<
"# looppoints for cell " << cellI <<
nl;
2111 pointField pointsOfLoop = loopPoints(loop, loopWeights);
2122 str <<
' ' << i + 1;
2124 str <<
' ' << 1 <<
nl;
2127 bool okLoop =
false;
2129 if (validEdgeLoop(loop, loopWeights))
2132 Map<edge> faceSplitCuts(loop.size());
2137 validLoop(cellI, loop, loopWeights, faceSplitCuts, anchorPoints);
2142 cellLoops_[cellI] = loop;
2143 cellAnchorPoints_[cellI].transfer(anchorPoints);
2148 faceSplitCut_.insert(iter.key(), iter());
2155 label cut = loop[cutI];
2159 label edgeI = getEdge(cut);
2161 edgeIsCut_[edgeI] =
true;
2163 edgeWeight_[edgeI] = loopWeights[cutI];
2167 label vertI = getVertex(cut);
2169 pointIsCut_[vertI] =
true;
2181 void Foam::cellCuts::setFromCellLoops
2185 const List<scalarField>& cellLoopWeights
2189 pointIsCut_ =
false;
2193 forAll(cellLabels, cellLabelI)
2195 label cellI = cellLabels[cellLabelI];
2197 const labelList& loop = cellLoops[cellLabelI];
2201 const scalarField& loopWeights = cellLoopWeights[cellLabelI];
2203 if (setFromCellLoop(cellI, loop, loopWeights))
2219 void Foam::cellCuts::setFromCellCutter
2221 const cellLooper& cellCutter,
2222 const List<refineCell>& refCells
2226 pointIsCut_ =
false;
2235 DynamicList<label> invalidCutCells(2);
2236 DynamicList<labelList> invalidCutLoops(2);
2237 DynamicList<scalarField> invalidCutLoopWeights(2);
2240 forAll(refCells, refCellI)
2242 const refineCell& refCell = refCells[refCellI];
2244 label cellI = refCell.cellNo();
2246 const vector& refDir = refCell.direction();
2268 if (setFromCellLoop(cellI, cellLoop, cellLoopWeights))
2274 cellLoops_[cellI].setSize(0);
2279 invalidCutCells.append(cellI);
2280 invalidCutLoops.append(cellLoop);
2281 invalidCutLoopWeights.append(cellLoopWeights);
2288 cellLoops_[cellI].setSize(0);
2292 if (debug && invalidCutCells.size())
2294 invalidCutCells.shrink();
2295 invalidCutLoops.shrink();
2296 invalidCutLoopWeights.shrink();
2298 fileName cutsFile(
"invalidLoopCells.obj");
2300 Pout<<
"cellCuts : writing inValidLoops cells to " << cutsFile <<
endl;
2302 OFstream cutsStream(cutsFile);
2313 fileName loopsFile(
"invalidLoops.obj");
2315 Pout<<
"cellCuts : writing inValidLoops loops to " << loopsFile <<
endl;
2317 OFstream loopsStream(loopsFile);
2321 forAll(invalidCutLoops, i)
2326 loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
2335 void Foam::cellCuts::setFromCellCutter
2337 const cellLooper& cellCutter,
2339 const List<plane>& cellCutPlanes
2343 pointIsCut_ =
false;
2352 DynamicList<label> invalidCutCells(2);
2353 DynamicList<labelList> invalidCutLoops(2);
2354 DynamicList<scalarField> invalidCutLoopWeights(2);
2359 label cellI = cellLabels[i];
2380 if (setFromCellLoop(cellI, cellLoop, cellLoopWeights))
2386 cellLoops_[cellI].setSize(0);
2391 invalidCutCells.append(cellI);
2392 invalidCutLoops.append(cellLoop);
2393 invalidCutLoopWeights.append(cellLoopWeights);
2400 cellLoops_[cellI].setSize(0);
2404 if (debug && invalidCutCells.size())
2406 invalidCutCells.shrink();
2407 invalidCutLoops.shrink();
2408 invalidCutLoopWeights.shrink();
2410 fileName cutsFile(
"invalidLoopCells.obj");
2412 Pout<<
"cellCuts : writing inValidLoops cells to " << cutsFile <<
endl;
2414 OFstream cutsStream(cutsFile);
2425 fileName loopsFile(
"invalidLoops.obj");
2427 Pout<<
"cellCuts : writing inValidLoops loops to " << loopsFile <<
endl;
2429 OFstream loopsStream(loopsFile);
2433 forAll(invalidCutLoops, i)
2438 loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
2447 void Foam::cellCuts::orientPlanesAndLoops()
2450 forAll(cellLoops_, cellI)
2452 const labelList& loop = cellLoops_[cellI];
2454 if (loop.size() && cellAnchorPoints_[cellI].empty())
2462 cellAnchorPoints_[cellI]
2469 Pout<<
"cellAnchorPoints:" <<
endl;
2471 forAll(cellAnchorPoints_, cellI)
2473 if (cellLoops_[cellI].size())
2475 if (cellAnchorPoints_[cellI].empty())
2478 <<
"No anchor points for cut cell " << cellI << endl
2484 Pout<<
" cell:" << cellI <<
" anchored at "
2485 << cellAnchorPoints_[cellI] <<
endl;
2493 forAll(cellLoops_, cellI)
2495 if (cellLoops_[cellI].size())
2504 void Foam::cellCuts::calcLoopsAndAddressing(
const labelList& cutCells)
2507 forAll(edgeIsCut_, edgeI)
2509 if (edgeIsCut_[edgeI])
2511 scalar weight = edgeWeight_[edgeI];
2513 if (weight < 0 || weight > 1)
2517 "cellCuts::calcLoopsAndAddressing(const labelList&)"
2518 ) <<
"Weight out of range [0,1]. Edge " << edgeI
2526 edgeWeight_[edgeI] = -GREAT;
2532 calcCellLoops(cutCells);
2537 forAll(cellLoops_, cellI)
2539 const labelList& loop = cellLoops_[cellI];
2543 Pout<<
"cell:" << cellI <<
" ";
2544 writeCuts(
Pout, loop, loopWeights(loop));
2556 void Foam::cellCuts::check()
const
2559 forAll(edgeIsCut_, edgeI)
2561 if (edgeIsCut_[edgeI])
2572 <<
"edge:" << edgeI <<
" vertices:"
2574 <<
" weight:" << edgeWeight_[edgeI] <<
" should have been"
2575 <<
" snapped to one of its endpoints"
2582 if (edgeWeight_[edgeI] > - 1)
2585 <<
"edge:" << edgeI <<
" vertices:"
2587 <<
" weight:" << edgeWeight_[edgeI] <<
" is not cut but"
2588 <<
" its weight is not marked invalid"
2595 forAll(cellLoops_, cellI)
2597 const labelList& loop = cellLoops_[cellI];
2601 label cut = loop[i];
2605 (isEdge(cut) && !edgeIsCut_[getEdge(cut)])
2606 || (!isEdge(cut) && !pointIsCut_[getVertex(cut)])
2610 writeCuts(
Pout, cuts, loopWeights(cuts));
2613 <<
"cell:" << cellI <<
" loop:"
2615 <<
" cut:" << cut <<
" is not marked as cut"
2622 forAll(cellLoops_, cellI)
2624 const labelList& anchors = cellAnchorPoints_[cellI];
2626 const labelList& loop = cellLoops_[cellI];
2628 if (loop.size() && anchors.empty())
2631 <<
"cell:" << cellI <<
" loop:" << loop
2632 <<
" has no anchor points"
2639 label cut = loop[i];
2644 &&
findIndex(anchors, getVertex(cut)) != -1
2648 <<
"cell:" << cellI <<
" loop:" << loop
2649 <<
" anchor points:" << anchors
2650 <<
" anchor:" << getVertex(cut) <<
" is part of loop"
2660 label faceI = iter.key();
2662 if (
mesh().isInternalFace(faceI))
2667 if (cellLoops_[own].empty() && cellLoops_[nei].empty())
2670 <<
"Internal face:" << faceI <<
" cut by " << iter()
2671 <<
" has owner:" << own
2672 <<
" and neighbour:" << nei
2673 <<
" that are both uncut"
2681 if (cellLoops_[own].empty())
2684 <<
"Boundary face:" << faceI <<
" cut by " << iter()
2685 <<
" has owner:" << own
2697 Foam::cellCuts::cellCuts
2707 pointIsCut_(expand(mesh.
nPoints(), meshVerts)),
2708 edgeIsCut_(expand(mesh.
nEdges(), meshEdges)),
2709 edgeWeight_(expand(mesh.
nEdges(), meshEdges, meshEdgeWeights)),
2711 faceSplitCut_(cutCells.
size()),
2712 cellLoops_(mesh.
nCells()),
2714 cellAnchorPoints_(mesh.
nCells())
2718 Pout<<
"cellCuts : constructor from cut verts and edges" <<
endl;
2721 calcLoopsAndAddressing(cutCells);
2724 orientPlanesAndLoops();
2735 Pout<<
"cellCuts : leaving constructor from cut verts and edges"
2743 Foam::cellCuts::cellCuts
2752 pointIsCut_(expand(mesh.
nPoints(), meshVerts)),
2753 edgeIsCut_(expand(mesh.
nEdges(), meshEdges)),
2754 edgeWeight_(expand(mesh.
nEdges(), meshEdges, meshEdgeWeights)),
2756 faceSplitCut_(mesh.
nFaces()/10 + 1),
2757 cellLoops_(mesh.
nCells()),
2759 cellAnchorPoints_(mesh.
nCells())
2763 Pout<<
"cellCuts : constructor from cellLoops" <<
endl;
2769 orientPlanesAndLoops();
2780 Pout<<
"cellCuts : leaving constructor from cellLoops" <<
endl;
2787 Foam::cellCuts::cellCuts
2796 pointIsCut_(mesh.
nPoints(),
false),
2797 edgeIsCut_(mesh.
nEdges(),
false),
2798 edgeWeight_(mesh.
nEdges(), -GREAT),
2800 faceSplitCut_(cellLabels.
size()),
2801 cellLoops_(mesh.
nCells()),
2803 cellAnchorPoints_(mesh.
nCells())
2807 Pout<<
"cellCuts : constructor from cellLoops" <<
endl;
2812 setFromCellLoops(cellLabels, cellLoops, cellEdgeWeights);
2815 orientPlanesAndLoops();
2826 Pout<<
"cellCuts : leaving constructor from cellLoops" <<
endl;
2832 Foam::cellCuts::cellCuts
2840 pointIsCut_(mesh.
nPoints(),
false),
2841 edgeIsCut_(mesh.
nEdges(),
false),
2842 edgeWeight_(mesh.
nEdges(), -GREAT),
2844 faceSplitCut_(refCells.
size()),
2845 cellLoops_(mesh.
nCells()),
2847 cellAnchorPoints_(mesh.
nCells())
2851 Pout<<
"cellCuts : constructor from cellCutter" <<
endl;
2856 setFromCellCutter(cellCutter, refCells);
2859 orientPlanesAndLoops();
2870 Pout<<
"cellCuts : leaving constructor from cellCutter" <<
endl;
2876 Foam::cellCuts::cellCuts
2885 pointIsCut_(mesh.
nPoints(),
false),
2886 edgeIsCut_(mesh.
nEdges(),
false),
2887 edgeWeight_(mesh.
nEdges(), -GREAT),
2889 faceSplitCut_(cellLabels.
size()),
2890 cellLoops_(mesh.
nCells()),
2892 cellAnchorPoints_(mesh.
nCells())
2896 Pout<<
"cellCuts : constructor from cellCutter with prescribed plane"
2902 setFromCellCutter(cellCutter, cellLabels, cutPlanes);
2905 orientPlanesAndLoops();
2916 Pout<<
"cellCuts : leaving constructor from cellCutter with prescribed"
2917 <<
" plane" <<
endl;
2923 Foam::cellCuts::cellCuts
2936 pointIsCut_(pointIsCut),
2937 edgeIsCut_(edgeIsCut),
2938 edgeWeight_(edgeWeight),
2940 faceSplitCut_(faceSplitCut),
2941 cellLoops_(cellLoops),
2943 cellAnchorPoints_(cellAnchorPoints)
2947 Pout<<
"cellCuts : constructor from components" <<
endl;
2948 Pout<<
"cellCuts : leaving constructor from components" <<
endl;
2971 const labelList& loop = cellLoops_[cellI];
2977 label cut = loop[fp];
2981 label edgeI = getEdge(cut);
2983 loopPts[fp] = coord(cut, edgeWeight_[edgeI]);
2987 loopPts[fp] = coord(cut, -GREAT);
3002 cellAnchorPoints_[cellI] =
3005 mesh().cellPoints()[cellI],
3006 cellAnchorPoints_[cellI],
3028 label startVertI = vertI;
3032 const point& pt = loopPts[fp];
3034 os <<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z() <<
endl;
3042 os <<
' ' << startVertI + fp + 1;
3052 forAll(cellLoops_, cellI)
3054 writeOBJ(os, loopPoints(cellI), vertI);
3061 const labelList& anchors = cellAnchorPoints_[cellI];
3063 writeOBJ(dir, cellI, loopPoints(cellI), anchors);