61 scalar faceAreaLimit = SMALL;
65 if (
mag(fAreas[facei]) > faceAreaLimit)
67 checkFaces.append(facei);
73 motionSmoother::checkMesh
101 for (
const label facei : badFaces)
108 syncTools::syncPointList
137 const label edgeI = faceEdges[eI];
138 const edge&
e = edges[edgeI];
158 edgeLabels[
count++] = edgeI;
162 if (
count != edgeLabels.size())
164 edgeLabels.setSize(
count);
171 void Foam::edgeCollapser::collapseToEdge
176 const vector& collapseAxis,
184 Map<point>& collapsePointToLocation
193 DynamicList<label> maxPriorityPts(
max(dNeg.size(), dPos.size()));
197 const label facePointi = facePtsNeg[fPtI];
198 const label facePtPriority = pointPriority[facePointi];
200 if (facePtPriority > maxPriority)
202 maxPriority = facePtPriority;
203 maxPriorityPts.clear();
204 maxPriorityPts.append(facePointi);
206 else if (facePtPriority == maxPriority)
208 maxPriorityPts.append(facePointi);
212 if (!maxPriorityPts.empty())
216 forAll(maxPriorityPts, ptI)
218 averagePt += pts[maxPriorityPts[ptI]];
221 collapseToPtA = averagePt/maxPriorityPts.size();
226 maxPriorityPts.clear();
228 labelList faceEdgesNeg = edgesFromPoints(facei, facePtsNeg);
234 collapsePointToLocation.set(facePtsNeg[pI], collapseToPtA);
244 const label facePointi = facePtsPos[fPtI];
245 const label facePtPriority = pointPriority[facePointi];
247 if (facePtPriority > maxPriority)
249 maxPriority = facePtPriority;
250 maxPriorityPts.clear();
251 maxPriorityPts.append(facePointi);
253 else if (facePtPriority == maxPriority)
255 maxPriorityPts.append(facePointi);
259 if (!maxPriorityPts.empty())
263 forAll(maxPriorityPts, ptI)
265 averagePt += pts[maxPriorityPts[ptI]];
268 collapseToPtB = averagePt/maxPriorityPts.size();
272 labelList faceEdgesPos = edgesFromPoints(facei, facePtsPos);
278 collapsePointToLocation.set(facePtsPos[pI], collapseToPtB);
283 void Foam::edgeCollapser::collapseToPoint
291 Map<point>& collapsePointToLocation
294 const face&
f = mesh_.faces()[facei];
299 DynamicList<label> maxPriorityPts(
f.size());
303 const label facePointi = facePts[fPtI];
304 const label facePtPriority = pointPriority[facePointi];
306 if (facePtPriority > maxPriority)
308 maxPriority = facePtPriority;
309 maxPriorityPts.clear();
310 maxPriorityPts.append(facePointi);
312 else if (facePtPriority == maxPriority)
314 maxPriorityPts.append(facePointi);
318 if (!maxPriorityPts.empty())
322 forAll(maxPriorityPts, ptI)
324 averagePt += pts[maxPriorityPts[ptI]];
327 collapseToPt = averagePt/maxPriorityPts.
size();
370 const labelList& faceEdges = mesh_.faceEdges()[facei];
376 collapsePointToLocation.set(
f[pI], collapseToPt);
381 void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio
398 scalar magJ =
mag(J);
413 detJ =
max(
det(J), SMALL);
420 collapseAxis =
f.edge(
f.longestEdge(pts)).unitVec(pts);
430 if (
mag(eVals.y() - eVals.x()) < 100*SMALL)
437 collapseAxis =
f.edge(
f.longestEdge(pts)).unitVec(pts);
465 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
467 const labelList& cellOwner = mesh_.faceOwner();
468 const labelList& cellNeighbour = mesh_.faceNeighbour();
470 const label nBoundaryFaces = mesh_.nBoundaryFaces();
473 for (label intFacei = 0; intFacei < mesh_.nInternalFaces(); ++intFacei)
475 const scalar cellOwnerVol =
max(0.0, V[cellOwner[intFacei]]);
476 const scalar cellNeighbourVol =
max(0.0, V[cellNeighbour[intFacei]]);
478 scalar targetFaceSizeA =
Foam::cbrt(cellOwnerVol);
479 scalar targetFaceSizeB =
Foam::cbrt(cellNeighbourVol);
481 targetFaceSizes[intFacei] = 0.5*(targetFaceSizeA + targetFaceSizeB);
491 label bFacei =
patch.start() - mesh_.nInternalFaces();
501 neiCellVolumes[bFacei++] =
max(0.0, V[faceCells[facei]]);
510 const label extFacei = patchFacei +
patch.start();
511 const scalar cellOwnerVol =
max(0.0, V[cellOwner[extFacei]]);
513 targetFaceSizes[extFacei] =
Foam::cbrt(cellOwnerVol);
524 label bFacei =
patch.start() - mesh_.nInternalFaces();
530 const label localFacei = patchFacei +
patch.start();
531 const scalar cellOwnerVol =
max(0.0, V[cellOwner[localFacei]]);
532 const scalar cellNeighbourVol = neiCellVolumes[bFacei++];
534 scalar targetFaceSizeA =
Foam::cbrt(cellOwnerVol);
535 scalar targetFaceSizeB =
Foam::cbrt(cellNeighbourVol);
537 targetFaceSizes[localFacei]
538 = 0.5*(targetFaceSizeA + targetFaceSizeB);
544 return targetFaceSizes;
553 const scalar targetFaceSize,
555 Map<point>& collapsePointToLocation,
559 const scalar collapseSizeLimitCoeff = faceFilterFactor[facei];
567 const scalar fA =
f.mag(pts);
570 scalar aspectRatio = 1.0;
572 faceCollapseAxisAndAspectRatio(
f, fC, collapseAxis, aspectRatio);
583 d[fPtI] = (collapseAxis & (pt - fC));
591 oldToNew =
invert(oldToNew.size(), oldToNew);
600 scalar dShift = -0.5*(d.first() + d.last());
630 SubList<scalar> dNeg(d, middle, 0);
631 SubList<label> facePtsNeg(facePts, middle, 0);
634 SubList<scalar> dPos(d, d.size() - middle, middle);
635 SubList<label> facePtsPos(facePts, d.size() - middle, middle);
647 if (dNeg.size() == 0 || dPos.size() == 0)
650 <<
"All points on one side of face centre, not collapsing."
658 collapseType typeOfCollapse = noCollapse;
660 if (
magSqr(collapseAxis) < VSMALL)
664 else if (fA < aspectRatio*
sqr(targetFaceSize*collapseSizeLimitCoeff))
668 allowEarlyCollapseToPoint_
669 && (d.last() - d.first())
671 *allowEarlyCollapseCoeff_*maxCollapseFaceToPointSideLengthCoeff_
678 (dNeg.last() < guardFraction_*dNeg.first())
679 && (dPos.first() > guardFraction_*dPos.last())
682 typeOfCollapse = toEdge;
686 (d.last() - d.first())
688 *maxCollapseFaceToPointSideLengthCoeff_
707 collapsePointToLocation
710 else if (typeOfCollapse == toEdge)
725 collapsePointToLocation
729 return typeOfCollapse;
733 Foam::label Foam::edgeCollapser::edgeMaster
739 label masterPoint = -1;
741 const label e0 =
e.start();
742 const label e1 =
e.end();
744 const label e0Priority = pointPriority[e0];
745 const label e1Priority = pointPriority[e1];
747 if (e0Priority > e1Priority)
751 else if (e0Priority < e1Priority)
755 else if (e0Priority == e1Priority)
791 void Foam::edgeCollapser::checkBoundaryPointMergeEdges
794 const label otherPointi,
796 Map<point>& collapsePointToLocation
801 const label e0Priority = pointPriority[pointi];
802 const label e1Priority = pointPriority[otherPointi];
804 if (e0Priority > e1Priority)
806 collapsePointToLocation.set
812 else if (e0Priority < e1Priority)
814 collapsePointToLocation.set
822 collapsePointToLocation.set
839 Foam::label Foam::edgeCollapser::breakStringsAtEdges
841 const bitSet& markedEdges,
843 List<pointEdgeCollapse>& allPointInfo
846 const edgeList& edges = mesh_.edges();
849 label nUncollapsed = 0;
853 if (markedEdges.test(eI))
855 const edge&
e = edges[eI];
857 const label startCollapseIndex
858 = allPointInfo[
e.start()].collapseIndex();
860 if (startCollapseIndex != -1 && startCollapseIndex != -2)
862 const label endCollapseIndex
863 = allPointInfo[
e.end()].collapseIndex();
868 && startCollapseIndex == endCollapseIndex
871 const labelList& ptEdgesStart = pointEdges[
e.start()];
873 forAll(ptEdgesStart, ptEdgeI)
875 const label edgeI = ptEdgesStart[ptEdgeI];
877 const label nbrPointi
878 = edges[edgeI].otherVertex(
e.start());
880 = allPointInfo[nbrPointi].collapseIndex();
884 nbrIndex == startCollapseIndex
901 void Foam::edgeCollapser::determineDuplicatePointsOnFace
904 bitSet& markedPoints,
907 List<pointEdgeCollapse>& allPointInfo
910 uniqueCollapses.clear();
911 duplicateCollapses.clear();
915 label index = allPointInfo[
f[fpI]].collapseIndex();
918 if (index != allPointInfo[
f.prevLabel(fpI)].collapseIndex())
920 if (!uniqueCollapses.insert(index))
923 duplicateCollapses.insert(index);
932 label index = allPointInfo[
f[fpI]].collapseIndex();
933 if (duplicateCollapses.found(index))
935 markedPoints.set(
f[fpI]);
941 Foam::label Foam::edgeCollapser::countEdgesOnFace
944 List<pointEdgeCollapse>& allPointInfo
951 const label pointi =
f[fpI];
952 const label
newPointi = allPointInfo[pointi].collapseIndex();
960 const label prevPointi =
f[
f.fcIndex(fpI)];
961 const label prevNewPointi
962 = allPointInfo[prevPointi].collapseIndex();
975 bool Foam::edgeCollapser::isFaceCollapsed
978 List<pointEdgeCollapse>& allPointInfo
981 label nEdges = countEdgesOnFace(
f, allPointInfo);
999 Foam::label Foam::edgeCollapser::syncCollapse
1001 const globalIndex& globalPoints,
1004 const Map<point>& collapsePointToLocation,
1005 List<pointEdgeCollapse>& allPointInfo
1008 const edgeList& edges = mesh_.edges();
1010 label nCollapsed = 0;
1012 DynamicList<label> initPoints(mesh_.nPoints());
1013 DynamicList<pointEdgeCollapse> initPointInfo(mesh_.nPoints());
1015 allPointInfo.
clear();
1016 allPointInfo.setSize(mesh_.nPoints());
1019 List<pointEdgeCollapse> allEdgeInfo
1022 pointEdgeCollapse(
Zero, -1, -1)
1030 const edge&
e = edges[edgeI];
1032 label masterPointi =
e.start();
1035 if (pointPriority[
e.end()] > pointPriority[
e.start()])
1037 masterPointi =
e.end();
1040 label masterPointPriority = pointPriority[masterPointi];
1042 label index = globalPoints.toGlobal(masterPointi);
1044 if (!collapsePointToLocation.found(masterPointi))
1046 const label otherVertex =
e.otherVertex(masterPointi);
1048 if (!collapsePointToLocation.found(otherVertex))
1051 << masterPointi <<
" on edge " << edgeI <<
" " <<
e
1052 <<
" is not marked for collapse."
1057 masterPointi = otherVertex;
1058 masterPointPriority = pointPriority[masterPointi];
1059 index = globalPoints.toGlobal(masterPointi);
1063 const point& collapsePoint = collapsePointToLocation[masterPointi];
1065 const pointEdgeCollapse pec
1074 allEdgeInfo[edgeI] = pointEdgeCollapse
1081 initPointInfo.append(pec);
1082 initPoints.append(
e.start());
1084 initPointInfo.append(pec);
1085 initPoints.append(
e.end());
1091 PointEdgeWave<pointEdgeCollapse> collapsePropagator
1098 mesh_.globalData().nTotalPoints()
1105 void Foam::edgeCollapser::filterFace
1107 const Map<DynamicList<label>>& collapseStrings,
1108 const List<pointEdgeCollapse>& allPointInfo,
1118 label pointi =
f[fp];
1120 label collapseIndex = allPointInfo[pointi].collapseIndex();
1123 if (collapseStrings.found(collapseIndex))
1125 const label localPointi = collapseStrings[collapseIndex][0];
1127 if (!SubList<label>(
f, newFp).found(localPointi))
1129 f[newFp++] = localPointi;
1132 else if (collapseIndex == -1)
1135 <<
"Point " << pointi <<
" was not visited by PointEdgeWave"
1140 f[newFp++] = pointi;
1151 const label size = newFp;
1155 for (label fp = 2; fp < size; fp++)
1160 label pointi =
f[fp];
1163 const label index = SubList<label>(
f, fp).find(pointi);
1168 <<
"Removing consecutive duplicate vertex in face "
1172 else if (index == fp2)
1175 <<
"Removing non-consecutive duplicate vertex in face "
1180 else if (index != -1)
1183 <<
"Pinched face " <<
f <<
endl;
1184 f[newFp++] = pointi;
1188 f[newFp++] = pointi;
1202 maxCollapseFaceToPointSideLengthCoeff_(0),
1203 allowEarlyCollapseToPoint_(false),
1204 allowEarlyCollapseCoeff_(0)
1208 Foam::edgeCollapser::edgeCollapser
1219 maxCollapseFaceToPointSideLengthCoeff_
1223 allowEarlyCollapseToPoint_
1227 allowEarlyCollapseCoeff_
1234 Info<<
"Edge Collapser Settings:" <<
nl
1235 <<
" Guard Fraction = " << guardFraction_ <<
nl
1236 <<
" Max collapse face to point side length = "
1237 << maxCollapseFaceToPointSideLengthCoeff_ <<
nl
1238 <<
" " << (allowEarlyCollapseToPoint_ ?
"Allow" :
"Do not allow")
1239 <<
" early collapse to point" <<
nl
1240 <<
" Early collapse coeff = " << allowEarlyCollapseCoeff_
1255 const labelList& faceOwner = mesh_.faceOwner();
1256 const labelList& faceNeighbour = mesh_.faceNeighbour();
1298 bool meshChanged =
false;
1300 bitSet removedPoints(mesh_.nPoints());
1308 forAll(allPointInfo, pointi)
1310 label collapseIndex = allPointInfo[pointi].collapseIndex();
1312 if (collapseIndex != -1 && collapseIndex != -2)
1314 ++(nPerIndex(collapseIndex, 0));
1319 collapseStrings.
resize(2*nPerIndex.size());
1326 forAll(allPointInfo, pointi)
1328 const label collapseIndex = allPointInfo[pointi].collapseIndex();
1330 if (collapseIndex != -1 && collapseIndex != -2)
1332 collapseStrings[collapseIndex].
append(pointi);
1364 bitSet cellRemoved(mesh_.nCells(),
false);
1366 label nUnvisited = 0;
1367 label nUncollapsed = 0;
1368 label nCollapsed = 0;
1388 label
nPoints = allPointInfo.size();
1397 <<
indent <<
"Not visited : " << nUnvisited <<
nl
1398 <<
indent <<
"Not collapsed : " << nUncollapsed <<
nl
1399 <<
indent <<
"Collapsed : " << nCollapsed <<
nl
1407 filterFace(collapseStrings, allPointInfo, newFaces[facei]);
1411 label nCellCollapsed = 0;
1415 if (!cellRemoved.test(celli))
1419 label nFaces = cFaces.size();
1423 label facei = cFaces[i];
1425 if (newFaces[facei].size() < 3)
1431 Pout<<
"Cell:" << celli
1432 <<
" uses faces:" << cFaces
1433 <<
" of which too many are marked for removal:"
1440 if (newFaces[cFaces[j]].size() < 3)
1442 Pout<<
' '<< cFaces[j];
1447 cellRemoved.set(celli);
1462 Info<<
indent <<
"Collapsing " << nCellCollapsed <<
" cells" <<
endl;
1464 if (nCellCollapsed == 0)
1473 bitSet doneFace(mesh_.nFaces(),
false);
1477 bitSet usedPoint(mesh_.nPoints(),
false);
1479 forAll(cellRemoved, celli)
1481 if (cellRemoved.test(celli))
1490 const face&
f = newFaces[facei];
1498 doneFace.set(facei);
1508 forAll(usedPoint, pointi)
1510 if (!usedPoint.test(pointi))
1512 removedPoints.set(pointi);
1520 forAll(allPointInfo, pointi)
1522 const label collapseIndex = allPointInfo[pointi].collapseIndex();
1523 const point& collapsePoint = allPointInfo[pointi].collapsePoint();
1527 !removedPoints.test(pointi)
1528 && collapseIndex != -1
1529 && collapseIndex != -2
1547 forAll(allPointInfo, pointi)
1549 if (removedPoints.test(pointi))
1551 const labelList& changedFaces = pointFaces[pointi];
1553 forAll(changedFaces, changedFacei)
1555 label facei = changedFaces[changedFacei];
1557 if (doneFace.set(facei))
1564 bool zoneFlip =
false;
1574 label own = faceOwner[facei];
1578 if (mesh_.isInternalFace(facei))
1580 nei = faceNeighbour[facei];
1616 const bool allowCellCollapse
1621 const faceList faces = mesh_.faces();
1622 const edgeList& edges = mesh_.edges();
1632 label nUncollapsed = 0;
1650 collapsePointToLocation,
1656 bitSet isCollapsedFace(mesh_.nFaces());
1657 bitSet markedPoints(mesh_.nPoints());
1661 const face&
f = faces[facei];
1663 isCollapsedFace[facei] = isFaceCollapsed(
f, allPointInfo);
1665 if (!isCollapsedFace.test(facei))
1667 determineDuplicatePointsOnFace
1688 forAll(markedPoints, pointi)
1690 if (markedPoints.test(pointi))
1692 const label index = allPointInfo[pointi].collapseIndex();
1694 const labelList& ptEdges = pointEdges[pointi];
1698 const label edgeI = ptEdges[ptEdgeI];
1699 const label nbrPointi = edges[edgeI].otherVertex(pointi);
1700 const label nbrIndex
1701 = allPointInfo[nbrPointi].collapseIndex();
1712 bitSet markedEdges(mesh_.nEdges());
1714 if (!allowCellCollapse)
1721 label nFaces = cFaces.size();
1725 label facei = cFaces[fI];
1727 if (isCollapsedFace.test(facei))
1737 label facei = cFaces[fI];
1739 const labelList& fEdges = faceEdges[facei];
1744 label edgeI = fEdges[fEdgeI];
1751 markedEdges.
set(edgeI);
1755 isCollapsedFace.unset(facei);
1763 <<
"Cell " << celli <<
" " << cFaces <<
nl
1764 <<
"is " << nFaces <<
", "
1765 <<
"but cell collapse has been disabled."
1779 nUncollapsed += breakStringsAtEdges
1788 Info<<
" Uncollapsed edges = " << nUncollapsed <<
" / "
1791 if (nUncollapsed == 0)
1811 const edgeList& edges = mesh_.edges();
1813 label nCollapsed = 0;
1817 const edge&
e = edges[edgeI];
1821 if (
e.mag(
points) < minEdgeLen[edgeI])
1825 label masterPointi = edgeMaster(pointPriority,
e);
1827 if (masterPointi == -1)
1832 collapsePointToLocation.set(
e.start(),
average);
1838 collapsePointToLocation.set(masterPointi, collapsePt);
1852 const scalar maxCos,
1858 const edgeList& edges = mesh_.edges();
1867 label nTotRemove = pointRemover.
countPointUsage(maxCos, pointCanBeDeleted);
1871 label nCollapsed = 0;
1875 forAll(pointEdges, pointi)
1877 if (pointCanBeDeleted[pointi])
1879 const labelList& pEdges = pointEdges[pointi];
1881 if (pEdges.size() == 2)
1885 label e0 = pEdges[0];
1886 label e1 = pEdges[1];
1891 scalar e0length =
mag
1896 scalar e1length =
mag
1901 if (e0length <= e1length)
1905 checkBoundaryPointMergeEdges
1908 edges[e0].otherVertex(pointi),
1910 collapsePointToLocation
1917 checkBoundaryPointMergeEdges
1920 edges[e1].otherVertex(pointi),
1922 collapsePointToLocation
1945 const faceList& faces = mesh_.faces();
1947 const scalarField targetFaceSizes = calcTargetFaceSizes();
1950 label nCollapseToPoint = 0;
1951 label nCollapseToEdge = 0;
1955 const face&
f = faces[fI];
1957 if (faceFilterFactor[fI] <= 0)
1967 targetFaceSizes[fI],
1969 collapsePointToLocation,
1973 if (flagCollapseFace == noCollapse)
1977 else if (flagCollapseFace ==
toPoint)
1981 else if (flagCollapseFace == toEdge)
1988 <<
"Face is marked to be collapsed to " << flagCollapseFace
1989 <<
". Currently can only collapse to point/edge."
1994 return labelPair(nCollapseToPoint, nCollapseToEdge);
2007 const faceList& faces = mesh_.faces();
2009 const scalarField targetFaceSizes = calcTargetFaceSizes();
2012 label nCollapseToPoint = 0;
2013 label nCollapseToEdge = 0;
2022 const face&
f = faces[fI];
2024 if (faceFilterFactor[fI] <= 0)
2034 targetFaceSizes[fI],
2036 collapsePointToLocation,
2040 if (flagCollapseFace == noCollapse)
2044 else if (flagCollapseFace ==
toPoint)
2048 else if (flagCollapseFace == toEdge)
2055 <<
"Face is marked to be collapsed to " << flagCollapseFace
2056 <<
". Currently can only collapse to point/edge."
2061 return labelPair(nCollapseToPoint, nCollapseToEdge);