59 x = (x==
y) ? x : value;
67 void Foam::hexRef8::reorder
88 newElems[newI] = elems[i];
92 elems.transfer(newElems);
96 void Foam::hexRef8::getFaceInfo
106 if (!mesh_.isInternalFace(faceI))
108 patchID = mesh_.boundaryMesh().whichPatch(faceI);
111 zoneID = mesh_.faceZones().whichZone(faceI);
117 const faceZone& fZone = mesh_.faceZones()[zoneID];
119 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
125 Foam::label Foam::hexRef8::addFace
127 polyTopoChange& meshMod,
134 label patchID, zoneID, zoneFlip;
136 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
140 if ((nei == -1) || (own < nei))
143 newFaceI = meshMod.setAction
163 newFaceI = meshMod.setAction
167 newFace.reverseFace(),
190 Foam::label Foam::hexRef8::addInternalFace
192 polyTopoChange& meshMod,
193 const label meshFaceI,
194 const label meshPointI,
200 if (mesh_.isInternalFace(meshFaceI))
202 return meshMod.setAction
230 return meshMod.setAction
286 void Foam::hexRef8::modFace
288 polyTopoChange& meshMod,
295 label patchID, zoneID, zoneFlip;
297 getFaceInfo(faceI, patchID, zoneID, zoneFlip);
301 (own != mesh_.faceOwner()[faceI])
303 mesh_.isInternalFace(faceI)
304 && (nei != mesh_.faceNeighbour()[faceI])
306 || (newFace != mesh_.faces()[faceI])
309 if ((nei == -1) || (own < nei))
333 newFace.reverseFace(),
350 Foam::scalar Foam::hexRef8::getLevel0EdgeLength()
const
352 if (cellLevel_.size() != mesh_.nCells())
356 "hexRef8::getLevel0EdgeLength() const"
357 ) <<
"Number of cells in mesh:" << mesh_.nCells()
358 <<
" does not equal size of cellLevel:" << cellLevel_.size()
360 <<
"This might be because of a restart with inconsistent cellLevel."
367 const scalar GREAT2 =
sqr(GREAT);
369 label nLevels =
gMax(cellLevel_)+1;
384 const label cLevel = cellLevel_[cellI];
386 const labelList& cEdges = mesh_.cellEdges(cellI);
390 label edgeI = cEdges[i];
392 if (edgeLevel[edgeI] == -1)
394 edgeLevel[edgeI] = cLevel;
396 else if (edgeLevel[edgeI] == labelMax)
400 else if (edgeLevel[edgeI] != cLevel)
402 edgeLevel[edgeI] = labelMax;
414 ifEqEqOp<labelMax>(),
423 const label eLevel = edgeLevel[edgeI];
425 if (eLevel >= 0 && eLevel < labelMax)
427 const edge&
e = mesh_.edges()[edgeI];
429 scalar edgeLenSqr =
magSqr(e.vec(mesh_.points()));
431 typEdgeLenSqr[eLevel] =
min(typEdgeLenSqr[eLevel], edgeLenSqr);
443 Pout<<
"hexRef8::getLevel0EdgeLength() :"
444 <<
" After phase1: Edgelengths (squared) per refinementlevel:"
445 << typEdgeLenSqr <<
endl;
459 const label cLevel = cellLevel_[cellI];
461 const labelList& cEdges = mesh_.cellEdges(cellI);
465 const edge&
e = mesh_.edges()[cEdges[i]];
467 scalar edgeLenSqr =
magSqr(e.vec(mesh_.points()));
469 maxEdgeLenSqr[cLevel] =
max(maxEdgeLenSqr[cLevel], edgeLenSqr);
478 Pout<<
"hexRef8::getLevel0EdgeLength() :"
479 <<
" Crappy Edgelengths (squared) per refinementlevel:"
480 << maxEdgeLenSqr <<
endl;
487 forAll(typEdgeLenSqr, levelI)
489 if (typEdgeLenSqr[levelI] == GREAT2 && maxEdgeLenSqr[levelI] >= 0)
491 typEdgeLenSqr[levelI] = maxEdgeLenSqr[levelI];
497 Pout<<
"hexRef8::getLevel0EdgeLength() :"
498 <<
" Final Edgelengths (squared) per refinementlevel:"
499 << typEdgeLenSqr <<
endl;
503 scalar level0Size = -1;
505 forAll(typEdgeLenSqr, levelI)
507 scalar lenSqr = typEdgeLenSqr[levelI];
515 Pout<<
"hexRef8::getLevel0EdgeLength() :"
516 <<
" For level:" << levelI
518 <<
" with equivalent level0 len:" << level0Size
525 if (level0Size == -1)
537 Foam::label Foam::hexRef8::getAnchorCell
546 if (cellAnchorPoints[cellI].size())
548 label index =
findIndex(cellAnchorPoints[cellI], pointI);
552 return cellAddedCells[cellI][index];
559 const face&
f = mesh_.faces()[faceI];
563 label index =
findIndex(cellAnchorPoints[cellI], f[fp]);
567 return cellAddedCells[cellI][index];
573 Perr<<
"cell:" << cellI <<
" anchorPoints:" << cellAnchorPoints[cellI]
577 <<
"Could not find point " << pointI
578 <<
" in the anchorPoints for cell " << cellI << endl
579 <<
"Does your original mesh obey the 2:1 constraint and"
580 <<
" did you use consistentRefinement to make your cells to refine"
581 <<
" obey this constraint as well?"
594 void Foam::hexRef8::getFaceNeighbours
610 mesh_.faceOwner()[faceI],
615 if (mesh_.isInternalFace(faceI))
621 mesh_.faceNeighbour()[faceI],
634 Foam::label Foam::hexRef8::findMinLevel(
const labelList& f)
const
636 label minLevel = labelMax;
641 label level = pointLevel_[f[fp]];
643 if (level < minLevel)
655 Foam::label Foam::hexRef8::findMaxLevel(
const labelList& f)
const
657 label maxLevel = labelMin;
662 label level = pointLevel_[f[fp]];
664 if (level > maxLevel)
675 Foam::label Foam::hexRef8::countAnchors
678 const label anchorLevel
685 if (pointLevel_[f[fp]] <= anchorLevel)
694 void Foam::hexRef8::dumpCell(
const label cellI)
const
696 OFstream str(mesh_.time().path()/
"cell_" +
Foam::name(cellI) +
".obj");
697 Pout<<
"hexRef8 : Dumping cell as obj to " << str.
name() <<
endl;
699 const cell& cFaces = mesh_.cells()[cellI];
701 Map<label> pointToObjVert;
706 const face& f = mesh_.faces()[cFaces[i]];
710 if (pointToObjVert.insert(f[fp], objVertI))
720 const face& f = mesh_.faces()[cFaces[i]];
724 label pointI = f[fp];
725 label nexPointI = f[f.fcIndex(fp)];
727 str <<
"l " << pointToObjVert[pointI]+1
728 <<
' ' << pointToObjVert[nexPointI]+1 <<
nl;
735 Foam::label Foam::hexRef8::findLevel
740 const bool searchForward,
741 const label wantedLevel
748 label pointI = f[fp];
750 if (pointLevel_[pointI] < wantedLevel)
752 dumpCell(mesh_.faceOwner()[faceI]);
753 if (mesh_.isInternalFace(faceI))
755 dumpCell(mesh_.faceNeighbour()[faceI]);
760 <<
" level:" << UIndirectList<label>(pointLevel_,
f)()
761 <<
" startFp:" << startFp
762 <<
" wantedLevel:" << wantedLevel
765 else if (pointLevel_[pointI] == wantedLevel)
780 dumpCell(mesh_.faceOwner()[faceI]);
781 if (mesh_.isInternalFace(faceI))
783 dumpCell(mesh_.faceNeighbour()[faceI]);
788 <<
" level:" << UIndirectList<label>(pointLevel_,
f)()
789 <<
" startFp:" << startFp
790 <<
" wantedLevel:" << wantedLevel
800 const face& f = mesh_.faces()[faceI];
804 return pointLevel_[f[findMaxLevel(f)]];
808 label ownLevel = cellLevel_[mesh_.faceOwner()[faceI]];
810 if (countAnchors(f, ownLevel) == 4)
814 else if (countAnchors(f, ownLevel+1) == 4)
826 void Foam::hexRef8::checkInternalOrientation
839 vector n(compactFace.normal(compactPoints));
841 vector dir(neiPt - ownPt);
846 <<
"cell:" << cellI <<
" old face:" << faceI
847 <<
" newFace:" << newFace <<
endl
848 <<
" coords:" << compactPoints
849 <<
" ownPt:" << ownPt
850 <<
" neiPt:" << neiPt
854 vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
856 scalar s = (fcToOwn&n) / (dir&n);
858 if (s < 0.1 || s > 0.9)
861 <<
"cell:" << cellI <<
" old face:" << faceI
862 <<
" newFace:" << newFace <<
endl
863 <<
" coords:" << compactPoints
864 <<
" ownPt:" << ownPt
865 <<
" neiPt:" << neiPt
872 void Foam::hexRef8::checkBoundaryOrientation
874 polyTopoChange& meshMod,
878 const point& boundaryPt,
882 face compactFace(
identity(newFace.size()));
883 pointField compactPoints(meshMod.points(), newFace);
885 vector n(compactFace.normal(compactPoints));
887 vector dir(boundaryPt - ownPt);
892 <<
"cell:" << cellI <<
" old face:" << faceI
893 <<
" newFace:" << newFace
894 <<
" coords:" << compactPoints
895 <<
" ownPt:" << ownPt
896 <<
" boundaryPt:" << boundaryPt
900 vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
902 scalar s = (fcToOwn&dir) /
magSqr(dir);
904 if (s < 0.7 || s > 1.3)
906 WarningIn(
"checkBoundaryOrientation(..)")
907 <<
"cell:" << cellI <<
" old face:" << faceI
908 <<
" newFace:" << newFace
909 <<
" coords:" << compactPoints
910 <<
" ownPt:" << ownPt
911 <<
" boundaryPt:" << boundaryPt
921 void Foam::hexRef8::insertEdgeSplit
926 DynamicList<label>& verts
929 if (p0 < mesh_.nPoints() && p1 < mesh_.nPoints())
933 if (edgeI != -1 && edgeMidPoint[edgeI] != -1)
935 verts.append(edgeMidPoint[edgeI]);
949 Foam::label Foam::hexRef8::storeMidPointInfo
957 const bool faceOrder,
958 const label edgeMidPointI,
959 const label anchorPointI,
960 const label faceMidPointI,
962 Map<edge>& midPointToAnchors,
963 Map<edge>& midPointToFaceMids,
964 polyTopoChange& meshMod
969 bool changed =
false;
970 bool haveTwoAnchors =
false;
974 if (edgeMidFnd == midPointToAnchors.end())
976 midPointToAnchors.insert(edgeMidPointI, edge(anchorPointI, -1));
980 edge&
e = edgeMidFnd();
982 if (anchorPointI != e[0])
991 if (e[0] != -1 && e[1] != -1)
993 haveTwoAnchors =
true;
997 bool haveTwoFaceMids =
false;
1001 if (faceMidFnd == midPointToFaceMids.end())
1003 midPointToFaceMids.insert(edgeMidPointI, edge(faceMidPointI, -1));
1007 edge& e = faceMidFnd();
1009 if (faceMidPointI != e[0])
1013 e[1] = faceMidPointI;
1018 if (e[0] != -1 && e[1] != -1)
1020 haveTwoFaceMids =
true;
1027 if (changed && haveTwoAnchors && haveTwoFaceMids)
1029 const edge& anchors = midPointToAnchors[edgeMidPointI];
1030 const edge& faceMids = midPointToFaceMids[edgeMidPointI];
1032 label otherFaceMidPointI = faceMids.otherVertex(faceMidPointI);
1039 DynamicList<label> newFaceVerts(4);
1040 if (faceOrder == (mesh_.faceOwner()[faceI] == cellI))
1042 newFaceVerts.append(faceMidPointI);
1053 newFaceVerts.append(edgeMidPointI);
1063 newFaceVerts.append(otherFaceMidPointI);
1064 newFaceVerts.append(cellMidPoint[cellI]);
1068 newFaceVerts.append(otherFaceMidPointI);
1078 newFaceVerts.append(edgeMidPointI);
1088 newFaceVerts.append(faceMidPointI);
1089 newFaceVerts.append(cellMidPoint[cellI]);
1093 newFace.transfer(newFaceVerts);
1095 label anchorCell0 = getAnchorCell
1103 label anchorCell1 = getAnchorCell
1109 anchors.otherVertex(anchorPointI)
1116 if (anchorCell0 < anchorCell1)
1121 ownPt = mesh_.points()[anchorPointI];
1122 neiPt = mesh_.points()[anchors.otherVertex(anchorPointI)];
1129 newFace = newFace.reverseFace();
1131 ownPt = mesh_.points()[anchors.otherVertex(anchorPointI)];
1132 neiPt = mesh_.points()[anchorPointI];
1139 if (anchorCell0 < anchorCell1)
1141 ownPt = mesh_.points()[anchorPointI];
1142 neiPt = mesh_.points()[anchors.otherVertex(anchorPointI)];
1146 ownPt = mesh_.points()[anchors.otherVertex(anchorPointI)];
1147 neiPt = mesh_.points()[anchorPointI];
1150 checkInternalOrientation
1161 return addInternalFace
1179 void Foam::hexRef8::createInternalFaces
1189 polyTopoChange& meshMod
1195 const cell& cFaces = mesh_.cells()[cellI];
1196 const label cLevel = cellLevel_[cellI];
1199 Map<edge> midPointToAnchors(24);
1201 Map<edge> midPointToFaceMids(24);
1204 DynamicList<label> storage;
1208 label nFacesAdded = 0;
1212 label faceI = cFaces[i];
1214 const face& f = mesh_.faces()[faceI];
1215 const labelList& fEdges = mesh_.faceEdges(faceI, storage);
1220 label faceMidPointI = -1;
1222 label nAnchors = countAnchors(f, cLevel);
1230 label anchorFp = -1;
1234 if (pointLevel_[f[fp]] <= cLevel)
1242 label edgeMid = findLevel
1246 f.fcIndex(anchorFp),
1250 label faceMid = findLevel
1259 faceMidPointI = f[faceMid];
1261 else if (nAnchors == 4)
1266 faceMidPointI = faceMidPoint[faceI];
1270 dumpCell(mesh_.faceOwner()[faceI]);
1271 if (mesh_.isInternalFace(faceI))
1273 dumpCell(mesh_.faceNeighbour()[faceI]);
1277 <<
"nAnchors:" << nAnchors
1278 <<
" faceI:" << faceI
1290 label point0 = f[fp0];
1292 if (pointLevel_[point0] <= cLevel)
1301 label edgeMidPointI = -1;
1303 label fp1 = f.fcIndex(fp0);
1305 if (pointLevel_[f[fp1]] <= cLevel)
1308 label edgeI = fEdges[fp0];
1310 edgeMidPointI = edgeMidPoint[edgeI];
1312 if (edgeMidPointI == -1)
1316 const labelList& cPoints = mesh_.cellPoints(cellI);
1319 <<
"cell:" << cellI <<
" cLevel:" << cLevel
1320 <<
" cell points:" << cPoints
1322 << UIndirectList<label>(pointLevel_, cPoints)()
1323 <<
" face:" << faceI
1326 << UIndirectList<label>(pointLevel_,
f)()
1327 <<
" faceAnchorLevel:" << faceAnchorLevel[faceI]
1328 <<
" faceMidPoint:" << faceMidPoint[faceI]
1329 <<
" faceMidPointI:" << faceMidPointI
1337 label edgeMid = findLevel(faceI, f, fp1,
true, cLevel+1);
1339 edgeMidPointI = f[edgeMid];
1342 label newFaceI = storeMidPointInfo
1365 if (nFacesAdded == 12)
1376 label fpMin1 = f.rcIndex(fp0);
1378 if (pointLevel_[f[fpMin1]] <= cLevel)
1381 label edgeI = fEdges[fpMin1];
1383 edgeMidPointI = edgeMidPoint[edgeI];
1385 if (edgeMidPointI == -1)
1389 const labelList& cPoints = mesh_.cellPoints(cellI);
1392 <<
"cell:" << cellI <<
" cLevel:" << cLevel
1393 <<
" cell points:" << cPoints
1395 << UIndirectList<label>(pointLevel_, cPoints)()
1396 <<
" face:" << faceI
1399 << UIndirectList<label>(pointLevel_,
f)()
1400 <<
" faceAnchorLevel:" << faceAnchorLevel[faceI]
1401 <<
" faceMidPoint:" << faceMidPoint[faceI]
1402 <<
" faceMidPointI:" << faceMidPointI
1410 label edgeMid = findLevel
1419 edgeMidPointI = f[edgeMid];
1422 newFaceI = storeMidPointInfo
1445 if (nFacesAdded == 12)
1453 if (nFacesAdded == 12)
1461 void Foam::hexRef8::walkFaceToMid
1466 const label startFp,
1467 DynamicList<label>& faceVerts
1470 const face& f = mesh_.faces()[faceI];
1471 const labelList& fEdges = mesh_.faceEdges(faceI);
1480 if (edgeMidPoint[fEdges[fp]] >= 0)
1482 faceVerts.
append(edgeMidPoint[fEdges[fp]]);
1487 if (pointLevel_[f[fp]] <= cLevel)
1493 else if (pointLevel_[f[fp]] == cLevel+1)
1496 faceVerts.append(f[fp]);
1500 else if (pointLevel_[f[fp]] == cLevel+2)
1503 faceVerts.append(f[fp]);
1510 void Foam::hexRef8::walkFaceFromMid
1515 const label startFp,
1516 DynamicList<label>& faceVerts
1519 const face& f = mesh_.faces()[faceI];
1520 const labelList& fEdges = mesh_.faceEdges(faceI);
1522 label fp = f.
rcIndex(startFp);
1526 if (pointLevel_[f[fp]] <= cLevel)
1531 else if (pointLevel_[f[fp]] == cLevel+1)
1534 faceVerts.append(f[fp]);
1537 else if (pointLevel_[f[fp]] == cLevel+2)
1547 if (edgeMidPoint[fEdges[fp]] >= 0)
1549 faceVerts.append(edgeMidPoint[fEdges[fp]]);
1558 faceVerts.append(f[fp]);
1565 Foam::label Foam::hexRef8::faceConsistentRefinement
1574 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1576 label own = mesh_.faceOwner()[faceI];
1577 label ownLevel = cellLevel_[own] + refineCell.get(own);
1579 label nei = mesh_.faceNeighbour()[faceI];
1580 label neiLevel = cellLevel_[nei] + refineCell.get(nei);
1582 if (ownLevel > (neiLevel+1))
1586 refineCell.set(nei, 1);
1590 refineCell.set(own, 0);
1594 else if (neiLevel > (ownLevel+1))
1598 refineCell.set(own, 1);
1602 refineCell.set(nei, 0);
1611 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1615 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1617 neiLevel[i] = cellLevel_[own] + refineCell.get(own);
1626 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1627 label ownLevel = cellLevel_[own] + refineCell.get(own);
1629 if (ownLevel > (neiLevel[i]+1))
1633 refineCell.set(own, 0);
1637 else if (neiLevel[i] > (ownLevel+1))
1641 refineCell.set(own, 1);
1652 void Foam::hexRef8::checkWantedRefinementLevels
1660 refineCell.set(cellsToRefine[i], 1);
1663 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1665 label own = mesh_.faceOwner()[faceI];
1666 label ownLevel = cellLevel_[own] + refineCell.get(own);
1668 label nei = mesh_.faceNeighbour()[faceI];
1669 label neiLevel = cellLevel_[nei] + refineCell.get(nei);
1671 if (
mag(ownLevel-neiLevel) > 1)
1677 "hexRef8::checkWantedRefinementLevels(const labelList&)"
1679 <<
" current level:" << cellLevel_[own]
1680 <<
" level after refinement:" << ownLevel
1682 <<
"neighbour cell:" << nei
1683 <<
" current level:" << cellLevel_[nei]
1684 <<
" level after refinement:" << neiLevel
1686 <<
"which does not satisfy 2:1 constraints anymore."
1693 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1697 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1699 neiLevel[i] = cellLevel_[own] + refineCell.get(own);
1708 label faceI = i + mesh_.nInternalFaces();
1710 label own = mesh_.faceOwner()[faceI];
1711 label ownLevel = cellLevel_[own] + refineCell.get(own);
1713 if (
mag(ownLevel - neiLevel[i]) > 1)
1715 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
1720 "hexRef8::checkWantedRefinementLevels(const labelList&)"
1721 ) <<
"Celllevel does not satisfy 2:1 constraint."
1722 <<
" On coupled face "
1724 <<
" on patch " << patchI <<
" "
1725 << mesh_.boundaryMesh()[patchI].name()
1726 <<
" owner cell " << own
1727 <<
" current level:" << cellLevel_[own]
1728 <<
" level after refinement:" << ownLevel
1730 <<
" (coupled) neighbour cell will get refinement "
1743 Pout<<
"hexRef8::setInstance(const fileName& inst) : "
1744 <<
"Resetting file instance to " << inst <<
endl;
1747 cellLevel_.instance() = inst;
1748 pointLevel_.instance() = inst;
1749 history_.instance() = inst;
1764 mesh_.facesInstance(),
1777 mesh_.facesInstance(),
1785 level0Edge_(getLevel0EdgeLength()),
1790 "refinementHistory",
1791 mesh_.facesInstance(),
1799 faceRemover_(mesh_, GREAT),
1800 savedPointLevel_(0),
1807 "hexRef8::hexRef8(const polyMesh&)"
1808 ) <<
"History enabled but number of visible cells in it "
1810 <<
" is not equal to the number of cells in the mesh "
1823 "hexRef8::hexRef8(const polyMesh&)"
1824 ) <<
"Restarted from inconsistent cellLevel or pointLevel files."
1826 <<
"Number of cells in mesh:" << mesh_.
nCells()
1827 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl
1828 <<
"Number of points in mesh:" << mesh_.
nPoints()
1829 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
1848 Foam::hexRef8::hexRef8
1862 mesh_.facesInstance(),
1875 mesh_.facesInstance(),
1883 level0Edge_(getLevel0EdgeLength()),
1888 "refinementHistory",
1889 mesh_.facesInstance(),
1897 faceRemover_(mesh_, GREAT),
1898 savedPointLevel_(0),
1901 if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
1905 "hexRef8::hexRef8(const polyMesh&, const labelList&"
1906 ", const labelList&, const refinementHistory&)"
1907 ) <<
"History enabled but number of visible cells in it "
1908 << history_.visibleCells().size()
1909 <<
" is not equal to the number of cells in the mesh "
1915 cellLevel_.size() != mesh_.nCells()
1916 || pointLevel_.size() != mesh_.nPoints()
1921 "hexRef8::hexRef8(const polyMesh&, const labelList&"
1922 ", const labelList&, const refinementHistory&)"
1923 ) <<
"Incorrect cellLevel or pointLevel size." <<
endl
1924 <<
"Number of cells in mesh:" << mesh_.nCells()
1925 <<
" does not equal size of cellLevel:" << cellLevel_.size() <<
endl
1926 <<
"Number of points in mesh:" << mesh_.nPoints()
1927 <<
" does not equal size of pointLevel:" << pointLevel_.size()
1932 checkRefinementLevels(-1,
labelList(0));
1945 Foam::hexRef8::hexRef8
1958 mesh_.facesInstance(),
1971 mesh_.facesInstance(),
1979 level0Edge_(getLevel0EdgeLength()),
1984 "refinementHistory",
1985 mesh_.facesInstance(),
1994 faceRemover_(mesh_, GREAT),
1995 savedPointLevel_(0),
2000 cellLevel_.size() != mesh_.nCells()
2001 || pointLevel_.size() != mesh_.nPoints()
2006 "hexRef8::hexRef8(const polyMesh&, const labelList&"
2007 ", const labelList&)"
2008 ) <<
"Incorrect cellLevel or pointLevel size." <<
endl
2009 <<
"Number of cells in mesh:" << mesh_.nCells()
2010 <<
" does not equal size of cellLevel:" << cellLevel_.size() <<
endl
2011 <<
"Number of points in mesh:" << mesh_.nPoints()
2012 <<
" does not equal size of pointLevel:" << pointLevel_.size()
2017 checkRefinementLevels(-1,
labelList(0));
2045 refineCell.set(cellsToRefine[i], 1);
2050 label nChanged = faceConsistentRefinement(maxSet, refineCell);
2056 Pout<<
"hexRef8::consistentRefinement : Changed " << nChanged
2057 <<
" refinement levels due to 2:1 conflicts."
2071 forAll(refineCell, cellI)
2073 if (refineCell.get(cellI) == 1)
2082 forAll(refineCell, cellI)
2084 if (refineCell.get(cellI) == 1)
2086 newCellsToRefine[nRefined++] = cellI;
2092 checkWantedRefinementLevels(newCellsToRefine);
2095 return newCellsToRefine;
2107 const label maxFaceDiff,
2110 const label maxPointDiff,
2114 const labelList& faceOwner = mesh_.faceOwner();
2115 const labelList& faceNeighbour = mesh_.faceNeighbour();
2118 if (maxFaceDiff <= 0)
2122 "hexRef8::consistentSlowRefinement"
2123 "(const label, const labelList&, const labelList&"
2124 ", const label, const labelList&)"
2125 ) <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl
2144 forAll(allCellInfo, cellI)
2150 maxFaceDiff*(cellLevel_[cellI]+1),
2151 maxFaceDiff*cellLevel_[cellI]
2158 label cellI = cellsToRefine[i];
2160 allCellInfo[cellI].count() = allCellInfo[cellI].refinementCount();
2175 label faceI = facesToCheck[i];
2177 if (allFaceInfo[faceI].valid())
2182 "hexRef8::consistentSlowRefinement"
2183 "(const label, const labelList&, const labelList&"
2184 ", const label, const labelList&)"
2185 ) <<
"Argument facesToCheck seems to have duplicate entries!"
2187 <<
"face:" << faceI <<
" occurs at positions "
2195 if (mesh_.isInternalFace(faceI))
2200 const refinementData& neiData = allCellInfo[faceNeighbour[faceI]];
2203 label faceRefineCount;
2206 faceCount = neiData.
count() + maxFaceDiff;
2207 faceRefineCount = faceCount + maxFaceDiff;
2211 faceCount = ownData.
count() + maxFaceDiff;
2212 faceRefineCount = faceCount + maxFaceDiff;
2215 seedFaces.append(faceI);
2216 seedFacesInfo.append
2224 allFaceInfo[faceI] = seedFacesInfo[seedFacesInfo.size()-1];
2228 label faceCount = ownData.
count() + maxFaceDiff;
2229 label faceRefineCount = faceCount + maxFaceDiff;
2231 seedFaces.append(faceI);
2232 seedFacesInfo.append
2240 allFaceInfo[faceI] = seedFacesInfo[seedFacesInfo.size()-1];
2248 forAll(faceNeighbour, faceI)
2251 if (!allFaceInfo[faceI].valid())
2253 label own = faceOwner[faceI];
2254 label nei = faceNeighbour[faceI];
2258 if (allCellInfo[own].count() > allCellInfo[nei].count())
2260 allFaceInfo[faceI].updateFace
2269 seedFacesInfo.append(allFaceInfo[faceI]);
2271 else if (allCellInfo[own].count() < allCellInfo[nei].count())
2273 allFaceInfo[faceI].updateFace
2281 seedFaces.append(faceI);
2282 seedFacesInfo.append(allFaceInfo[faceI]);
2290 for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
2293 if (!allFaceInfo[faceI].valid())
2295 label own = faceOwner[faceI];
2307 seedFaces.append(faceI);
2308 seedFacesInfo.append(faceData);
2325 Pout<<
"hexRef8::consistentSlowRefinement : Seeded "
2326 << seedFaces.size() <<
" faces between cells with different"
2327 <<
" refinement level." <<
endl;
2331 levelCalc.
setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
2333 seedFacesInfo.clear();
2336 levelCalc.
iterate(mesh_.globalData().nTotalFaces());
2344 if (maxPointDiff == -1)
2352 labelList maxPointCount(mesh_.nPoints(), 0);
2354 forAll(maxPointCount, pointI)
2356 label& pLevel = maxPointCount[pointI];
2358 const labelList& pCells = mesh_.pointCells(pointI);
2362 pLevel =
max(pLevel, allCellInfo[pCells[i]].count());
2384 label pointI = pointsToCheck[i];
2389 const labelList& pCells = mesh_.pointCells(pointI);
2393 label cellI = pCells[pCellI];
2401 maxPointCount[pointI]
2402 > cellInfo.
count() + maxFaceDiff*maxPointDiff
2410 const cell& cFaces = mesh_.cells()[cellI];
2414 label faceI = cFaces[cFaceI];
2426 if (faceData.
count() > allFaceInfo[faceI].count())
2428 changedFacesInfo.insert(faceI, faceData);
2435 label nChanged = changedFacesInfo.size();
2445 seedFaces.setCapacity(changedFacesInfo.size());
2446 seedFacesInfo.setCapacity(changedFacesInfo.size());
2450 seedFaces.append(iter.key());
2451 seedFacesInfo.append(iter());
2458 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
2460 label own = mesh_.faceOwner()[faceI];
2463 + (allCellInfo[own].isRefined() ? 1 : 0);
2465 label nei = mesh_.faceNeighbour()[faceI];
2468 + (allCellInfo[nei].isRefined() ? 1 : 0);
2470 if (
mag(ownLevel-neiLevel) > 1)
2476 "hexRef8::consistentSlowRefinement"
2477 "(const label, const labelList&, const labelList&"
2478 ", const label, const labelList&)"
2480 <<
" current level:" << cellLevel_[own]
2481 <<
" current refData:" << allCellInfo[own]
2482 <<
" level after refinement:" << ownLevel
2484 <<
"neighbour cell:" << nei
2485 <<
" current level:" << cellLevel_[nei]
2486 <<
" current refData:" << allCellInfo[nei]
2487 <<
" level after refinement:" << neiLevel
2489 <<
"which does not satisfy 2:1 constraints anymore." <<
nl
2490 <<
"face:" << faceI <<
" faceRefData:" << allFaceInfo[faceI]
2499 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2500 labelList neiCount(mesh_.nFaces()-mesh_.nInternalFaces());
2501 labelList neiRefCount(mesh_.nFaces()-mesh_.nInternalFaces());
2505 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
2506 neiLevel[i] = cellLevel_[own];
2507 neiCount[i] = allCellInfo[own].count();
2508 neiRefCount[i] = allCellInfo[own].refinementCount();
2519 label faceI = i+mesh_.nInternalFaces();
2521 label own = mesh_.faceOwner()[faceI];
2524 + (allCellInfo[own].isRefined() ? 1 : 0);
2528 + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
2530 if (
mag(ownLevel - nbrLevel) > 1)
2533 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
2537 "hexRef8::consistentSlowRefinement"
2538 "(const label, const labelList&, const labelList&"
2539 ", const label, const labelList&)"
2540 ) <<
"Celllevel does not satisfy 2:1 constraint."
2541 <<
" On coupled face "
2543 <<
" refData:" << allFaceInfo[faceI]
2544 <<
" on patch " << patchI <<
" "
2545 << mesh_.boundaryMesh()[patchI].name() <<
nl
2546 <<
"owner cell " << own
2547 <<
" current level:" << cellLevel_[own]
2548 <<
" current count:" << allCellInfo[own].count()
2549 <<
" current refCount:"
2550 << allCellInfo[own].refinementCount()
2551 <<
" level after refinement:" << ownLevel
2553 <<
"(coupled) neighbour cell"
2554 <<
" has current level:" << neiLevel[i]
2555 <<
" current count:" << neiCount[i]
2556 <<
" current refCount:" << neiRefCount[i]
2557 <<
" level after refinement:" << nbrLevel
2567 forAll(allCellInfo, cellI)
2569 if (allCellInfo[cellI].isRefined())
2579 forAll(allCellInfo, cellI)
2581 if (allCellInfo[cellI].isRefined())
2583 newCellsToRefine[nRefined++] = cellI;
2589 Pout<<
"hexRef8::consistentSlowRefinement : From "
2590 << cellsToRefine.
size() <<
" to " << newCellsToRefine.
size()
2591 <<
" cells to refine." <<
endl;
2594 return newCellsToRefine;
2600 const label maxFaceDiff,
2605 const labelList& faceOwner = mesh_.faceOwner();
2606 const labelList& faceNeighbour = mesh_.faceNeighbour();
2608 if (maxFaceDiff <= 0)
2612 "hexRef8::consistentSlowRefinement2"
2613 "(const label, const labelList&, const labelList&)"
2614 ) <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl
2618 const scalar level0Size = 2*maxFaceDiff*level0Edge_;
2636 label cellI = cellsToRefine[i];
2641 mesh_.cellCentres()[cellI],
2646 forAll(allCellInfo, cellI)
2648 if (!allCellInfo[cellI].valid())
2653 mesh_.cellCentres()[cellI],
2670 label faceI = facesToCheck[i];
2672 if (allFaceInfo[faceI].valid())
2677 "hexRef8::consistentSlowRefinement2"
2678 "(const label, const labelList&, const labelList&)"
2679 ) <<
"Argument facesToCheck seems to have duplicate entries!"
2681 <<
"face:" << faceI <<
" occurs at positions "
2686 label own = faceOwner[faceI];
2690 allCellInfo[own].valid()
2691 ? allCellInfo[own].originLevel()
2695 if (!mesh_.isInternalFace(faceI))
2699 const point& fc = mesh_.faceCentres()[faceI];
2708 allFaceInfo[faceI].updateFace
2719 label nei = faceNeighbour[faceI];
2723 allCellInfo[nei].valid()
2724 ? allCellInfo[nei].originLevel()
2728 if (ownLevel == neiLevel)
2731 allFaceInfo[faceI].updateFace
2739 allFaceInfo[faceI].updateFace
2751 allFaceInfo[faceI].updateFace
2759 allFaceInfo[faceI].updateFace
2770 seedFacesInfo.append(allFaceInfo[faceI]);
2777 forAll(faceNeighbour, faceI)
2780 if (!allFaceInfo[faceI].valid())
2782 label own = faceOwner[faceI];
2786 allCellInfo[own].valid()
2787 ? allCellInfo[own].originLevel()
2791 label nei = faceNeighbour[faceI];
2795 allCellInfo[nei].valid()
2796 ? allCellInfo[nei].originLevel()
2800 if (ownLevel > neiLevel)
2804 allFaceInfo[faceI].updateFace
2812 seedFacesInfo.append(allFaceInfo[faceI]);
2814 else if (neiLevel > ownLevel)
2816 seedFaces.append(faceI);
2817 allFaceInfo[faceI].updateFace
2825 seedFacesInfo.append(allFaceInfo[faceI]);
2831 seedFacesInfo.shrink();
2841 mesh_.globalData().nTotalCells()
2881 label cellI = cellsToRefine[i];
2883 allCellInfo[cellI].originLevel() =
sizeof(label)*8-2;
2884 allCellInfo[cellI].origin() = cc[cellI];
2891 forAll(allCellInfo, cellI)
2893 label wanted = allCellInfo[cellI].wantedLevel(cc[cellI]);
2895 if (wanted > cellLevel_[cellI]+1)
2897 refineCell.set(cellI, 1);
2900 faceConsistentRefinement(
true, refineCell);
2904 label nChanged = faceConsistentRefinement(
true, refineCell);
2910 Pout<<
"hexRef8::consistentSlowRefinement2 : Changed " << nChanged
2911 <<
" refinement levels due to 2:1 conflicts."
2924 forAll(refineCell, cellI)
2926 if (refineCell.get(cellI) == 1)
2935 forAll(refineCell, cellI)
2937 if (refineCell.get(cellI) == 1)
2939 newCellsToRefine[nRefined++] = cellI;
2945 Pout<<
"hexRef8::consistentSlowRefinement2 : From "
2946 << cellsToRefine.
size() <<
" to " << newCellsToRefine.
size()
2947 <<
" cells to refine." <<
endl;
2952 cellSet cellsIn(mesh_,
"cellsToRefineIn", cellsToRefine);
2953 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
2954 << cellsIn.
size() <<
" to cellSet "
2959 cellSet cellsOut(mesh_,
"cellsToRefineOut", newCellsToRefine);
2960 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
2961 << cellsOut.
size() <<
" to cellSet "
2968 forAll(newCellsToRefine, i)
2970 refineCell.set(newCellsToRefine[i], 1);
2974 label nChanged = faceConsistentRefinement(
true, refineCell);
2979 mesh_,
"cellsToRefineOut2", newCellsToRefine.
size()
2981 forAll(refineCell, cellI)
2983 if (refineCell.get(cellI) == 1)
2985 cellsOut2.insert(cellI);
2988 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
2989 << cellsOut2.size() <<
" to cellSet "
2990 << cellsOut2.objectPath() <<
endl;
2996 forAll(refineCell, cellI)
3000 refineCell.get(cellI) == 1
3001 && savedRefineCell.
get(cellI) == 0
3007 "hexRef8::consistentSlowRefinement2"
3008 "(const label, const labelList&, const labelList&)"
3009 ) <<
"Cell:" << cellI <<
" cc:"
3010 << mesh_.cellCentres()[cellI]
3011 <<
" was not marked for refinement but does not obey"
3012 <<
" 2:1 constraints."
3019 return newCellsToRefine;
3032 Pout<<
"hexRef8::setRefinement :"
3033 <<
" Checking initial mesh just to make sure" <<
endl;
3042 savedPointLevel_.clear();
3043 savedCellLevel_.clear();
3048 forAll(cellLevel_, cellI)
3050 newCellLevel.append(cellLevel_[cellI]);
3053 forAll(pointLevel_, pointI)
3055 newPointLevel.append(pointLevel_[pointI]);
3061 Pout<<
"hexRef8::setRefinement :"
3062 <<
" Allocating " << cellLabels.
size() <<
" cell midpoints."
3070 labelList cellMidPoint(mesh_.nCells(), -1);
3074 label cellI = cellLabels[i];
3076 label anchorPointI = mesh_.faces()[mesh_.cells()[cellI][0]][0];
3082 mesh_.cellCentres()[cellI],
3089 newPointLevel(cellMidPoint[cellI]) = cellLevel_[cellI]+1;
3095 cellSet splitCells(mesh_,
"splitCells", cellLabels.
size());
3097 forAll(cellMidPoint, cellI)
3099 if (cellMidPoint[cellI] >= 0)
3101 splitCells.insert(cellI);
3105 Pout<<
"hexRef8::setRefinement : Dumping " << splitCells.size()
3106 <<
" cells to split to cellSet " << splitCells.objectPath()
3119 Pout<<
"hexRef8::setRefinement :"
3120 <<
" Allocating edge midpoints."
3129 labelList edgeMidPoint(mesh_.nEdges(), -1);
3132 forAll(cellMidPoint, cellI)
3134 if (cellMidPoint[cellI] >= 0)
3136 const labelList& cEdges = mesh_.cellEdges(cellI);
3140 label edgeI = cEdges[i];
3142 const edge& e = mesh_.edges()[edgeI];
3146 pointLevel_[e[0]] <= cellLevel_[cellI]
3147 && pointLevel_[e[1]] <= cellLevel_[cellI]
3150 edgeMidPoint[edgeI] = 12345;
3178 forAll(edgeMidPoint, edgeI)
3180 if (edgeMidPoint[edgeI] >= 0)
3183 edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points());
3191 point(-GREAT, -GREAT, -GREAT),
3197 forAll(edgeMidPoint, edgeI)
3199 if (edgeMidPoint[edgeI] >= 0)
3204 const edge& e = mesh_.edges()[edgeI];
3217 newPointLevel(edgeMidPoint[edgeI]) =
3230 OFstream str(mesh_.time().path()/
"edgeMidPoint.obj");
3232 forAll(edgeMidPoint, edgeI)
3234 if (edgeMidPoint[edgeI] >= 0)
3236 const edge& e = mesh_.edges()[edgeI];
3242 Pout<<
"hexRef8::setRefinement :"
3243 <<
" Dumping edge centres to split to file " << str.
name() <<
endl;
3253 Pout<<
"hexRef8::setRefinement :"
3254 <<
" Allocating face midpoints."
3260 labelList faceAnchorLevel(mesh_.nFaces());
3262 for (label faceI = 0; faceI < mesh_.nFaces(); faceI++)
3264 faceAnchorLevel[faceI] = getAnchorLevel(faceI);
3269 labelList faceMidPoint(mesh_.nFaces(), -1);
3274 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
3276 if (faceAnchorLevel[faceI] >= 0)
3278 label own = mesh_.faceOwner()[faceI];
3279 label ownLevel = cellLevel_[own];
3280 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3282 label nei = mesh_.faceNeighbour()[faceI];
3283 label neiLevel = cellLevel_[nei];
3284 label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
3288 newOwnLevel > faceAnchorLevel[faceI]
3289 || newNeiLevel > faceAnchorLevel[faceI]
3292 faceMidPoint[faceI] = 12345;
3305 labelList newNeiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
3309 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
3310 label ownLevel = cellLevel_[own];
3311 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3313 newNeiLevel[i] = newOwnLevel;
3323 label faceI = i+mesh_.nInternalFaces();
3325 if (faceAnchorLevel[faceI] >= 0)
3327 label own = mesh_.faceOwner()[faceI];
3328 label ownLevel = cellLevel_[own];
3329 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3333 newOwnLevel > faceAnchorLevel[faceI]
3334 || newNeiLevel[i] > faceAnchorLevel[faceI]
3337 faceMidPoint[faceI] = 12345;
3363 mesh_.nFaces()-mesh_.nInternalFaces(),
3364 point(-GREAT, -GREAT, -GREAT)
3369 label faceI = i+mesh_.nInternalFaces();
3371 if (faceMidPoint[faceI] >= 0)
3373 bFaceMids[i] = mesh_.faceCentres()[faceI];
3384 forAll(faceMidPoint, faceI)
3386 if (faceMidPoint[faceI] >= 0)
3391 const face& f = mesh_.faces()[faceI];
3398 faceI < mesh_.nInternalFaces()
3399 ? mesh_.faceCentres()[faceI]
3400 : bFaceMids[faceI-mesh_.nInternalFaces()]
3410 newPointLevel(faceMidPoint[faceI]) = faceAnchorLevel[faceI]+1;
3417 faceSet splitFaces(mesh_,
"splitFaces", cellLabels.
size());
3419 forAll(faceMidPoint, faceI)
3421 if (faceMidPoint[faceI] >= 0)
3423 splitFaces.insert(faceI);
3427 Pout<<
"hexRef8::setRefinement : Dumping " << splitFaces.size()
3428 <<
" faces to split to faceSet " << splitFaces.objectPath() <<
endl;
3449 Pout<<
"hexRef8::setRefinement :"
3450 <<
" Finding cell anchorPoints (8 per cell)"
3461 labelList nAnchorPoints(mesh_.nCells(), 0);
3463 forAll(cellMidPoint, cellI)
3465 if (cellMidPoint[cellI] >= 0)
3467 cellAnchorPoints[cellI].setSize(8);
3471 forAll(pointLevel_, pointI)
3473 const labelList& pCells = mesh_.pointCells(pointI);
3477 label cellI = pCells[pCellI];
3481 cellMidPoint[cellI] >= 0
3482 && pointLevel_[pointI] <= cellLevel_[cellI]
3485 if (nAnchorPoints[cellI] == 8)
3490 "hexRef8::setRefinement(const labelList&"
3491 ", polyTopoChange&)"
3492 ) <<
"cell " << cellI
3493 <<
" of level " << cellLevel_[cellI]
3494 <<
" uses more than 8 points of equal or"
3495 <<
" lower level" <<
nl
3496 <<
"Points so far:" << cellAnchorPoints[cellI]
3500 cellAnchorPoints[cellI][nAnchorPoints[cellI]++] = pointI;
3505 forAll(cellMidPoint, cellI)
3507 if (cellMidPoint[cellI] >= 0)
3509 if (nAnchorPoints[cellI] != 8)
3513 const labelList& cPoints = mesh_.cellPoints(cellI);
3517 "hexRef8::setRefinement(const labelList&"
3518 ", polyTopoChange&)"
3519 ) <<
"cell " << cellI
3520 <<
" of level " << cellLevel_[cellI]
3521 <<
" does not seem to have 8 points of equal or"
3522 <<
" lower level" <<
endl
3523 <<
"cellPoints:" << cPoints <<
endl
3538 Pout<<
"hexRef8::setRefinement :"
3539 <<
" Adding cells (1 per anchorPoint)"
3546 forAll(cellAnchorPoints, cellI)
3548 const labelList& cAnchors = cellAnchorPoints[cellI];
3550 if (cAnchors.
size() == 8)
3552 labelList& cAdded = cellAddedCells[cellI];
3559 newCellLevel[cellI] = cellLevel_[cellI]+1;
3562 for (label i = 1; i < 8; i++)
3572 mesh_.cellZones().whichZone(cellI)
3576 newCellLevel(cAdded[i]) = cellLevel_[cellI]+1;
3591 Pout<<
"hexRef8::setRefinement :"
3592 <<
" Marking faces to be handled"
3600 forAll(cellMidPoint, cellI)
3602 if (cellMidPoint[cellI] >= 0)
3604 const cell& cFaces = mesh_.cells()[cellI];
3608 affectedFace.set(cFaces[i], 1);
3613 forAll(faceMidPoint, faceI)
3615 if (faceMidPoint[faceI] >= 0)
3617 affectedFace.set(faceI, 1);
3621 forAll(edgeMidPoint, edgeI)
3623 if (edgeMidPoint[edgeI] >= 0)
3625 const labelList& eFaces = mesh_.edgeFaces(edgeI);
3629 affectedFace.set(eFaces[i], 1);
3641 Pout<<
"hexRef8::setRefinement : Splitting faces" <<
endl;
3644 forAll(faceMidPoint, faceI)
3646 if (faceMidPoint[faceI] >= 0 && affectedFace.get(faceI) == 1)
3652 const face& f = mesh_.faces()[faceI];
3656 bool modifiedFace =
false;
3658 label anchorLevel = faceAnchorLevel[faceI];
3664 label pointI = f[fp];
3666 if (pointLevel_[pointI] <= anchorLevel)
3672 faceVerts.
append(pointI);
3688 faceVerts.
append(faceMidPoint[faceI]);
3721 if (mesh_.isInternalFace(faceI))
3723 label oldOwn = mesh_.faceOwner()[faceI];
3724 label oldNei = mesh_.faceNeighbour()[faceI];
3726 checkInternalOrientation
3731 mesh_.cellCentres()[oldOwn],
3732 mesh_.cellCentres()[oldNei],
3738 label oldOwn = mesh_.faceOwner()[faceI];
3740 checkBoundaryOrientation
3745 mesh_.cellCentres()[oldOwn],
3746 mesh_.faceCentres()[faceI],
3755 modifiedFace =
true;
3757 modFace(meshMod, faceI, newFace, own, nei);
3761 addFace(meshMod, faceI, newFace, own, nei);
3767 affectedFace.set(faceI, 0);
3777 Pout<<
"hexRef8::setRefinement :"
3778 <<
" Adding edge splits to unsplit faces"
3785 forAll(edgeMidPoint, edgeI)
3787 if (edgeMidPoint[edgeI] >= 0)
3791 const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
3795 label faceI = eFaces[i];
3797 if (faceMidPoint[faceI] < 0 && affectedFace.get(faceI) == 1)
3801 const face& f = mesh_.faces()[faceI];
3802 const labelList& fEdges = mesh_.faceEdges
3812 newFaceVerts.append(f[fp]);
3814 label edgeI = fEdges[fp];
3816 if (edgeMidPoint[edgeI] >= 0)
3818 newFaceVerts.
append(edgeMidPoint[edgeI]);
3827 label anchorFp = findMinLevel(f);
3844 if (mesh_.isInternalFace(faceI))
3846 label oldOwn = mesh_.faceOwner()[faceI];
3847 label oldNei = mesh_.faceNeighbour()[faceI];
3849 checkInternalOrientation
3854 mesh_.cellCentres()[oldOwn],
3855 mesh_.cellCentres()[oldNei],
3861 label oldOwn = mesh_.faceOwner()[faceI];
3863 checkBoundaryOrientation
3868 mesh_.cellCentres()[oldOwn],
3869 mesh_.faceCentres()[faceI],
3875 modFace(meshMod, faceI, newFace, own, nei);
3878 affectedFace.set(faceI, 0);
3890 Pout<<
"hexRef8::setRefinement :"
3891 <<
" Changing owner/neighbour for otherwise unaffected faces"
3895 forAll(affectedFace, faceI)
3897 if (affectedFace.get(faceI) == 1)
3899 const face& f = mesh_.faces()[faceI];
3903 label anchorFp = findMinLevel(f);
3917 modFace(meshMod, faceI, f, own, nei);
3920 affectedFace.set(faceI, 0);
3935 Pout<<
"hexRef8::setRefinement :"
3936 <<
" Create new internal faces for split cells"
3940 forAll(cellMidPoint, cellI)
3942 if (cellMidPoint[cellI] >= 0)
3964 label minPointI = labelMax;
3965 label maxPointI = labelMin;
3967 forAll(cellMidPoint, cellI)
3969 if (cellMidPoint[cellI] >= 0)
3971 minPointI =
min(minPointI, cellMidPoint[cellI]);
3972 maxPointI =
max(maxPointI, cellMidPoint[cellI]);
3975 forAll(faceMidPoint, faceI)
3977 if (faceMidPoint[faceI] >= 0)
3979 minPointI =
min(minPointI, faceMidPoint[faceI]);
3980 maxPointI =
max(maxPointI, faceMidPoint[faceI]);
3983 forAll(edgeMidPoint, edgeI)
3985 if (edgeMidPoint[edgeI] >= 0)
3987 minPointI =
min(minPointI, edgeMidPoint[edgeI]);
3988 maxPointI =
max(maxPointI, edgeMidPoint[edgeI]);
3992 if (minPointI != labelMax && minPointI != mesh_.nPoints())
3995 <<
"Added point labels not consecutive to existing mesh points."
3997 <<
"mesh_.nPoints():" << mesh_.nPoints()
3998 <<
" minPointI:" << minPointI
3999 <<
" maxPointI:" << maxPointI
4004 pointLevel_.transfer(newPointLevel);
4005 cellLevel_.transfer(newCellLevel);
4008 setInstance(mesh_.facesInstance());
4015 if (history_.active())
4019 Pout<<
"hexRef8::setRefinement :"
4020 <<
" Updating refinement history to " << cellLevel_.size()
4021 <<
" cells" <<
endl;
4025 history_.resize(cellLevel_.size());
4027 forAll(cellAddedCells, cellI)
4029 const labelList& addedCells = cellAddedCells[cellI];
4031 if (addedCells.
size())
4034 history_.storeSplit(cellI, addedCells);
4045 label cellI = cellLabels[i];
4047 refinedCells[i].
transfer(cellAddedCells[cellI]);
4050 return refinedCells;
4061 savedPointLevel_.resize(2*pointsToStore.
size());
4064 label pointI = pointsToStore[i];
4065 savedPointLevel_.insert(pointI, pointLevel_[pointI]);
4068 savedCellLevel_.
resize(2*cellsToStore.
size());
4071 label cellI = cellsToStore[i];
4072 savedCellLevel_.insert(cellI, cellLevel_[cellI]);
4084 updateMesh(map, dummyMap, dummyMap, dummyMap);
4102 Pout<<
"hexRef8::updateMesh :"
4103 <<
" Updating various lists"
4112 Pout<<
"hexRef8::updateMesh :"
4115 <<
" nCells:" << mesh_.nCells()
4117 <<
" cellLevel_:" << cellLevel_.size()
4120 <<
" nPoints:" << mesh_.nPoints()
4122 <<
" pointLevel_:" << pointLevel_.size()
4126 if (reverseCellMap.
size() == cellLevel_.size())
4132 reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
4140 forAll(cellMap, newCellI)
4142 label oldCellI = cellMap[newCellI];
4146 newCellLevel[newCellI] = -1;
4150 newCellLevel[newCellI] = cellLevel_[oldCellI];
4160 label newCellI = iter.key();
4161 label storedCellI = iter();
4165 if (fnd == savedCellLevel_.
end())
4167 FatalErrorIn(
"hexRef8::updateMesh(const mapPolyMesh&)")
4168 <<
"Problem : trying to restore old value for new cell "
4169 << newCellI <<
" but cannot find old cell " << storedCellI
4170 <<
" in map of stored values " << savedCellLevel_
4173 cellLevel_[newCellI] = fnd();
4194 if (reversePointMap.
size() == pointLevel_.size())
4197 reorder(reversePointMap, mesh_.nPoints(), -1, pointLevel_);
4206 forAll(pointMap, newPointI)
4208 label oldPointI = pointMap[newPointI];
4210 if (oldPointI == -1)
4219 newPointLevel[newPointI] = -1;
4223 newPointLevel[newPointI] = pointLevel_[oldPointI];
4226 pointLevel_.
transfer(newPointLevel);
4233 label newPointI = iter.key();
4234 label storedPointI = iter();
4238 if (fnd == savedPointLevel_.
end())
4240 FatalErrorIn(
"hexRef8::updateMesh(const mapPolyMesh&)")
4241 <<
"Problem : trying to restore old value for new point "
4242 << newPointI <<
" but cannot find old point "
4244 <<
" in map of stored values " << savedPointLevel_
4247 pointLevel_[newPointI] = fnd();
4264 if (history_.active())
4266 history_.updateMesh(map);
4270 setInstance(mesh_.facesInstance());
4273 faceRemover_.updateMesh(map);
4288 Pout<<
"hexRef8::subset :"
4289 <<
" Updating various lists"
4293 if (history_.active())
4297 "hexRef8::subset(const labelList&, const labelList&"
4298 ", const labelList&)"
4299 ) <<
"Subsetting will not work in combination with unrefinement."
4301 <<
"Proceed at your own risk." <<
endl;
4309 forAll(cellMap, newCellI)
4311 newCellLevel[newCellI] = cellLevel_[cellMap[newCellI]];
4314 cellLevel_.transfer(newCellLevel);
4320 <<
"cellLevel_ contains illegal value -1 after mapping:"
4330 forAll(pointMap, newPointI)
4332 newPointLevel[newPointI] = pointLevel_[pointMap[newPointI]];
4335 pointLevel_.transfer(newPointLevel);
4341 <<
"pointLevel_ contains illegal value -1 after mapping:"
4348 if (history_.active())
4350 history_.subset(pointMap, faceMap, cellMap);
4354 setInstance(mesh_.facesInstance());
4366 Pout<<
"hexRef8::distribute :"
4367 <<
" Updating various lists"
4377 if (history_.active())
4379 history_.distribute(map);
4383 faceRemover_.distribute(map);
4389 const scalar smallDim = 1
E-6 * mesh_.globalData().bb().mag();
4393 Pout<<
"hexRef8::checkMesh : Using matching tolerance smallDim:"
4394 << smallDim <<
endl;
4404 labelList nei(mesh_.nFaces()-mesh_.nInternalFaces());
4407 nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4426 label faceI = pp.
start();
4430 label own = mesh_.faceOwner()[faceI];
4431 label bFaceI = faceI-mesh_.nInternalFaces();
4437 <<
"Faces do not seem to be correct across coupled"
4438 <<
" boundaries" <<
endl
4439 <<
"Coupled face " << faceI
4440 <<
" between owner " << own
4441 <<
" on patch " << pp.
name()
4442 <<
" and coupled neighbour " << nei[bFaceI]
4443 <<
" has two faces connected to it:"
4459 scalarField neiFaceAreas(mesh_.nFaces()-mesh_.nInternalFaces());
4462 neiFaceAreas[i] =
mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
4470 label faceI = i+mesh_.nInternalFaces();
4472 const scalar magArea =
mag(mesh_.faceAreas()[faceI]);
4474 if (
mag(magArea - neiFaceAreas[i]) > smallDim)
4476 const face& f = mesh_.faces()[faceI];
4477 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4479 dumpCell(mesh_.faceOwner()[faceI]);
4482 <<
"Faces do not seem to be correct across coupled"
4483 <<
" boundaries" <<
endl
4484 <<
"Coupled face " << faceI
4485 <<
" on patch " << patchI
4486 <<
" " << mesh_.boundaryMesh()[patchI].
name()
4488 <<
" has face area:" << magArea
4489 <<
" (coupled) neighbour face area differs:"
4491 <<
" to within tolerance " << smallDim
4501 labelList nVerts(mesh_.nFaces()-mesh_.nInternalFaces());
4505 nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size();
4513 label faceI = i+mesh_.nInternalFaces();
4515 const face& f = mesh_.faces()[faceI];
4517 if (f.
size() != nVerts[i])
4519 dumpCell(mesh_.faceOwner()[faceI]);
4521 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4524 <<
"Faces do not seem to be correct across coupled"
4525 <<
" boundaries" <<
endl
4526 <<
"Coupled face " << faceI
4527 <<
" on patch " << patchI
4528 <<
" " << mesh_.boundaryMesh()[patchI].
name()
4530 <<
" has size:" << f.
size()
4531 <<
" (coupled) neighbour face has size:"
4543 pointField anchorPoints(mesh_.nFaces()-mesh_.nInternalFaces());
4547 label faceI = i+mesh_.nInternalFaces();
4548 const point& fc = mesh_.faceCentres()[faceI];
4549 const face& f = mesh_.faces()[faceI];
4550 const vector anchorVec(mesh_.points()[f[0]] - fc);
4552 anchorPoints[i] = anchorVec;
4562 label faceI = i+mesh_.nInternalFaces();
4563 const point& fc = mesh_.faceCentres()[faceI];
4564 const face& f = mesh_.faces()[faceI];
4565 const vector anchorVec(mesh_.points()[f[0]] - fc);
4567 if (
mag(anchorVec - anchorPoints[i]) > smallDim)
4569 dumpCell(mesh_.faceOwner()[faceI]);
4571 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4574 <<
"Faces do not seem to be correct across coupled"
4575 <<
" boundaries" <<
endl
4576 <<
"Coupled face " << faceI
4577 <<
" on patch " << patchI
4578 <<
" " << mesh_.boundaryMesh()[patchI].
name()
4580 <<
" has anchor vector:" << anchorVec
4581 <<
" (coupled) neighbour face anchor vector differs:"
4583 <<
" to within tolerance " << smallDim
4591 Pout<<
"hexRef8::checkMesh : Returning" <<
endl;
4598 const label maxPointDiff,
4604 Pout<<
"hexRef8::checkRefinementLevels :"
4605 <<
" Checking 2:1 refinement level" <<
endl;
4610 cellLevel_.size() != mesh_.nCells()
4611 || pointLevel_.size() != mesh_.nPoints()
4614 FatalErrorIn(
"hexRef8::checkRefinementLevels(const label)")
4615 <<
"cellLevel size should be number of cells"
4616 <<
" and pointLevel size should be number of points."<<
nl
4617 <<
"cellLevel:" << cellLevel_.size()
4618 <<
" mesh.nCells():" << mesh_.nCells() <<
nl
4619 <<
"pointLevel:" << pointLevel_.size()
4620 <<
" mesh.nPoints():" << mesh_.nPoints()
4630 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
4632 label own = mesh_.faceOwner()[faceI];
4633 label nei = mesh_.faceNeighbour()[faceI];
4635 if (
mag(cellLevel_[own] - cellLevel_[nei]) > 1)
4642 "hexRef8::checkRefinementLevels(const label)"
4643 ) <<
"Celllevel does not satisfy 2:1 constraint." <<
nl
4644 <<
"On face " << faceI <<
" owner cell " << own
4645 <<
" has refinement " << cellLevel_[own]
4646 <<
" neighbour cell " << nei <<
" has refinement "
4653 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
4657 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4659 neiLevel[i] = cellLevel_[own];
4667 label faceI = i+mesh_.nInternalFaces();
4669 label own = mesh_.faceOwner()[faceI];
4671 if (
mag(cellLevel_[own] - neiLevel[i]) > 1)
4675 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
4679 "hexRef8::checkRefinementLevels(const label)"
4680 ) <<
"Celllevel does not satisfy 2:1 constraint."
4681 <<
" On coupled face " << faceI
4682 <<
" on patch " << patchI <<
" "
4683 << mesh_.boundaryMesh()[patchI].name()
4684 <<
" owner cell " << own <<
" has refinement "
4686 <<
" (coupled) neighbour cell has refinement "
4710 forAll(syncPointLevel, pointI)
4712 if (pointLevel_[pointI] != syncPointLevel[pointI])
4716 "hexRef8::checkRefinementLevels(const label)"
4717 ) <<
"PointLevel is not consistent across coupled patches."
4719 <<
"point:" << pointI <<
" coord:" << mesh_.points()[pointI]
4720 <<
" has level " << pointLevel_[pointI]
4721 <<
" whereas the coupled point has level "
4722 << syncPointLevel[pointI]
4730 if (maxPointDiff != -1)
4733 labelList maxPointLevel(mesh_.nPoints(), 0);
4735 forAll(maxPointLevel, pointI)
4737 const labelList& pCells = mesh_.pointCells(pointI);
4739 label& pLevel = maxPointLevel[pointI];
4743 pLevel =
max(pLevel, cellLevel_[pCells[i]]);
4760 label pointI = pointsToCheck[i];
4762 const labelList& pCells = mesh_.pointCells(pointI);
4766 label cellI = pCells[i];
4770 mag(cellLevel_[cellI]-maxPointLevel[pointI])
4778 "hexRef8::checkRefinementLevels(const label)"
4779 ) <<
"Too big a difference between"
4780 <<
" point-connected cells." <<
nl
4782 <<
" cellLevel:" << cellLevel_[cellI]
4783 <<
" uses point:" << pointI
4784 <<
" coord:" << mesh_.points()[pointI]
4785 <<
" which is also used by a cell with level:"
4786 << maxPointLevel[pointI]
4872 checkRefinementLevels(-1,
labelList(0));
4877 Pout<<
"hexRef8::getSplitPoints :"
4878 <<
" Calculating unrefineable points" <<
endl;
4882 if (!history_.active())
4885 <<
"Only call if constructed with history capability"
4893 labelList splitMaster(mesh_.nPoints(), -1);
4894 labelList splitMasterLevel(mesh_.nPoints(), 0);
4899 for (label pointI = 0; pointI < mesh_.nPoints(); pointI++)
4901 const labelList& pCells = mesh_.pointCells(pointI);
4903 if (pCells.
size() != 8)
4905 splitMaster[pointI] = -2;
4910 const labelList& visibleCells = history_.visibleCells();
4912 forAll(visibleCells, cellI)
4914 const labelList& cPoints = mesh_.cellPoints(cellI);
4916 if (visibleCells[cellI] != -1 && history_.parentIndex(cellI) >= 0)
4918 label parentIndex = history_.parentIndex(cellI);
4923 label pointI = cPoints[i];
4925 label masterCellI = splitMaster[pointI];
4927 if (masterCellI == -1)
4934 splitMaster[pointI] = parentIndex;
4935 splitMasterLevel[pointI] = cellLevel_[cellI] - 1;
4937 else if (masterCellI == -2)
4943 (masterCellI != parentIndex)
4944 || (splitMasterLevel[pointI] != cellLevel_[cellI] - 1)
4949 splitMaster[pointI] = -2;
4958 label pointI = cPoints[i];
4960 splitMaster[pointI] = -2;
4968 label faceI = mesh_.nInternalFaces();
4969 faceI < mesh_.nFaces();
4973 const face& f = mesh_.faces()[faceI];
4977 splitMaster[f[fp]] = -2;
4984 label nSplitPoints = 0;
4986 forAll(splitMaster, pointI)
4988 if (splitMaster[pointI] >= 0)
4997 forAll(splitMaster, pointI)
4999 if (splitMaster[pointI] >= 0)
5001 splitPoints[nSplitPoints++] = pointI;
5079 Pout<<
"hexRef8::consistentUnrefinement :"
5080 <<
" Determining 2:1 consistent unrefinement" <<
endl;
5087 "hexRef8::consistentUnrefinement(const labelList&, const bool"
5088 ) <<
"maxSet not implemented yet."
5100 forAll(pointsToUnrefine, i)
5102 label pointI = pointsToUnrefine[i];
5104 unrefinePoint.set(pointI, 1);
5115 forAll(unrefinePoint, pointI)
5117 if (unrefinePoint.get(pointI) == 1)
5119 const labelList& pCells = mesh_.pointCells(pointI);
5123 unrefineCell.set(pCells[j], 1);
5136 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
5138 label own = mesh_.faceOwner()[faceI];
5139 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5141 label nei = mesh_.faceNeighbour()[faceI];
5142 label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
5144 if (ownLevel < (neiLevel-1))
5151 unrefineCell.set(nei, 1);
5155 if (unrefineCell.get(own) == 0)
5161 unrefineCell.set(own, 0);
5165 else if (neiLevel < (ownLevel-1))
5169 unrefineCell.set(own, 1);
5173 if (unrefineCell.get(nei) == 0)
5179 unrefineCell.set(nei, 0);
5187 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
5191 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
5193 neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
5201 label faceI = i+mesh_.nInternalFaces();
5202 label own = mesh_.faceOwner()[faceI];
5203 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5205 if (ownLevel < (neiLevel[i]-1))
5209 if (unrefineCell.get(own) == 0)
5215 unrefineCell.set(own, 0);
5219 else if (neiLevel[i] < (ownLevel-1))
5223 if (unrefineCell.get(own) == 1)
5229 unrefineCell.set(own, 1);
5239 Pout<<
"hexRef8::consistentUnrefinement :"
5240 <<
" Changed " << nChanged
5241 <<
" refinement levels due to 2:1 conflicts."
5255 forAll(unrefinePoint, pointI)
5257 if (unrefinePoint.get(pointI) == 1)
5259 const labelList& pCells = mesh_.pointCells(pointI);
5263 if (unrefineCell.get(pCells[j]) == 0)
5265 unrefinePoint.set(pointI, 0);
5277 forAll(unrefinePoint, pointI)
5279 if (unrefinePoint.get(pointI) == 1)
5288 forAll(unrefinePoint, pointI)
5290 if (unrefinePoint.get(pointI) == 1)
5292 newPointsToUnrefine[nSet++] = pointI;
5296 return newPointsToUnrefine;
5306 if (!history_.active())
5310 "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5311 ) <<
"Only call if constructed with history capability"
5317 Pout<<
"hexRef8::setUnrefinement :"
5318 <<
" Checking initial mesh just to make sure" <<
endl;
5322 forAll(cellLevel_, cellI)
5324 if (cellLevel_[cellI] < 0)
5328 "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5329 ) <<
"Illegal cell level " << cellLevel_[cellI]
5330 <<
" for cell " << cellI
5337 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5340 cellSet cSet(mesh_,
"splitPointCells", splitPointLabels.
size());
5342 forAll(splitPointLabels, i)
5344 const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
5348 cSet.insert(pCells[j]);
5353 Pout<<
"hexRef8::setRefinement : Dumping " << pSet.
size()
5355 << cSet.size() <<
" cells for unrefinement to" <<
nl
5357 <<
" cellSet " << cSet.objectPath()
5369 forAll(splitPointLabels, i)
5375 splitFaces.insert(pFaces[j]);
5381 faceRemover_.compatibleRemoves
5389 if (facesToRemove.
size() != splitFaces.size())
5393 "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5394 ) <<
"Ininitial set of split points to unrefine does not"
5395 <<
" seem to be consistent or not mid points of refined cells"
5404 forAll(splitPointLabels, i)
5406 label pointI = splitPointLabels[i];
5410 const labelList& pCells = mesh_.pointCells(pointI);
5413 if (pCells.
size() != 8)
5417 "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
5418 ) <<
"splitPoint " << pointI
5419 <<
" should have 8 cells using it. It has " << pCells
5428 label masterCellI =
min(pCells);
5432 label cellI = pCells[j];
5434 label region = cellRegion[cellI];
5439 <<
"Ininitial set of split points to unrefine does not"
5440 <<
" seem to be consistent or not mid points"
5441 <<
" of refined cells" <<
nl
5442 <<
"cell:" << cellI <<
" on splitPoint " << pointI
5443 <<
" has no region to be merged into"
5447 if (masterCellI != cellRegionMaster[region])
5450 <<
"cell:" << cellI <<
" on splitPoint:" << pointI
5451 <<
" in region " << region
5452 <<
" has master:" << cellRegionMaster[region]
5453 <<
" which is not the lowest numbered cell"
5454 <<
" among the pointCells:" << pCells
5463 faceRemover_.setRefinement
5474 forAll(splitPointLabels, i)
5476 label pointI = splitPointLabels[i];
5478 const labelList& pCells = mesh_.pointCells(pointI);
5480 label masterCellI =
min(pCells);
5484 cellLevel_[pCells[j]]--;
5487 history_.combineCells(masterCellI, pCells);
5491 setInstance(mesh_.facesInstance());
5500 bool writeOk = cellLevel_.write() && pointLevel_.write();
5502 if (history_.active())
5504 writeOk = writeOk && history_.write();