63 x = (
x ==
y) ?
x : value;
71 void Foam::hexRef8::reorder
92 newElems[newI] = elems[i];
96 elems.transfer(newElems);
100 void Foam::hexRef8::getFaceInfo
110 if (!mesh_.isInternalFace(facei))
112 patchID = mesh_.boundaryMesh().whichPatch(facei);
115 zoneID = mesh_.faceZones().whichZone(facei);
121 const faceZone& fZone = mesh_.faceZones()[
zoneID];
123 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
129 Foam::label Foam::hexRef8::addFace
131 polyTopoChange& meshMod,
144 if ((nei == -1) || (own < nei))
147 newFacei = meshMod.setAction
167 newFacei = meshMod.setAction
171 newFace.reverseFace(),
194 Foam::label Foam::hexRef8::addInternalFace
196 polyTopoChange& meshMod,
197 const label meshFacei,
198 const label meshPointi,
204 if (mesh_.isInternalFace(meshFacei))
206 return meshMod.setAction
234 return meshMod.setAction
290 void Foam::hexRef8::modFace
292 polyTopoChange& meshMod,
305 (own != mesh_.faceOwner()[facei])
307 mesh_.isInternalFace(facei)
308 && (nei != mesh_.faceNeighbour()[facei])
310 || (newFace != mesh_.faces()[facei])
313 if ((nei == -1) || (own < nei))
337 newFace.reverseFace(),
354 Foam::scalar Foam::hexRef8::getLevel0EdgeLength()
const
356 if (cellLevel_.size() != mesh_.nCells())
359 <<
"Number of cells in mesh:" << mesh_.nCells()
360 <<
" does not equal size of cellLevel:" << cellLevel_.size()
362 <<
"This might be because of a restart with inconsistent cellLevel."
369 const scalar GREAT2 =
sqr(GREAT);
371 label nLevels =
gMax(cellLevel_)+1;
386 const label cLevel = cellLevel_[celli];
388 const labelList& cEdges = mesh_.cellEdges(celli);
392 label edgeI = cEdges[i];
394 if (edgeLevel[edgeI] == -1)
396 edgeLevel[edgeI] = cLevel;
398 else if (edgeLevel[edgeI] ==
labelMax)
402 else if (edgeLevel[edgeI] != cLevel)
416 ifEqEqOp<labelMax>(),
424 const label eLevel = edgeLevel[edgeI];
426 if (eLevel >= 0 && eLevel <
labelMax)
428 const edge&
e = mesh_.edges()[edgeI];
430 scalar edgeLenSqr =
magSqr(
e.vec(mesh_.points()));
432 typEdgeLenSqr[eLevel] =
min(typEdgeLenSqr[eLevel], edgeLenSqr);
444 Pout<<
"hexRef8::getLevel0EdgeLength() :"
445 <<
" After phase1: Edgelengths (squared) per refinementlevel:"
446 << typEdgeLenSqr <<
endl;
460 const label cLevel = cellLevel_[celli];
462 const labelList& cEdges = mesh_.cellEdges(celli);
466 const edge&
e = mesh_.edges()[cEdges[i]];
468 scalar edgeLenSqr =
magSqr(
e.vec(mesh_.points()));
470 maxEdgeLenSqr[cLevel] =
max(maxEdgeLenSqr[cLevel], edgeLenSqr);
479 Pout<<
"hexRef8::getLevel0EdgeLength() :"
480 <<
" Poor Edgelengths (squared) per refinementlevel:"
481 << maxEdgeLenSqr <<
endl;
488 forAll(typEdgeLenSqr, levelI)
490 if (typEdgeLenSqr[levelI] == GREAT2 && maxEdgeLenSqr[levelI] >= 0)
492 typEdgeLenSqr[levelI] = maxEdgeLenSqr[levelI];
498 Pout<<
"hexRef8::getLevel0EdgeLength() :"
499 <<
" Final Edgelengths (squared) per refinementlevel:"
500 << typEdgeLenSqr <<
endl;
504 scalar level0Size = -1;
506 forAll(typEdgeLenSqr, levelI)
508 scalar lenSqr = typEdgeLenSqr[levelI];
516 Pout<<
"hexRef8::getLevel0EdgeLength() :"
517 <<
" For level:" << levelI
519 <<
" with equivalent level0 len:" << level0Size
526 if (level0Size == -1)
538 Foam::label Foam::hexRef8::getAnchorCell
547 if (cellAnchorPoints[celli].size())
549 label index = cellAnchorPoints[celli].find(pointi);
553 return cellAddedCells[celli][index];
560 const face&
f = mesh_.faces()[facei];
564 label index = cellAnchorPoints[celli].find(
f[fp]);
568 return cellAddedCells[celli][index];
574 Perr<<
"cell:" << celli <<
" anchorPoints:" << cellAnchorPoints[celli]
578 <<
"Could not find point " << pointi
579 <<
" in the anchorPoints for cell " << celli <<
endl
580 <<
"Does your original mesh obey the 2:1 constraint and"
581 <<
" did you use consistentRefinement to make your cells to refine"
582 <<
" obey this constraint as well?"
595 void Foam::hexRef8::getFaceNeighbours
611 mesh_.faceOwner()[facei],
616 if (mesh_.isInternalFace(facei))
622 mesh_.faceNeighbour()[facei],
635 Foam::label Foam::hexRef8::findMinLevel(
const labelList&
f)
const
642 label level = pointLevel_[
f[fp]];
644 if (level < minLevel)
656 Foam::label Foam::hexRef8::findMaxLevel(
const labelList&
f)
const
663 label level = pointLevel_[
f[fp]];
665 if (level > maxLevel)
676 Foam::label Foam::hexRef8::countAnchors
679 const label anchorLevel
686 if (pointLevel_[
f[fp]] <= anchorLevel)
695 void Foam::hexRef8::dumpCell(
const label celli)
const
697 OFstream str(mesh_.time().path()/
"cell_" +
Foam::name(celli) +
".obj");
698 Pout<<
"hexRef8 : Dumping cell as obj to " << str.
name() <<
endl;
700 const cell& cFaces = mesh_.cells()[celli];
702 Map<label> pointToObjVert;
707 const face&
f = mesh_.faces()[cFaces[i]];
711 if (pointToObjVert.insert(
f[fp], objVertI))
721 const face&
f = mesh_.faces()[cFaces[i]];
725 label pointi =
f[fp];
726 label nexPointi =
f[
f.fcIndex(fp)];
728 str <<
"l " << pointToObjVert[pointi]+1
729 <<
' ' << pointToObjVert[nexPointi]+1 <<
nl;
736 Foam::label Foam::hexRef8::findLevel
741 const bool searchForward,
742 const label wantedLevel
749 label pointi =
f[fp];
751 if (pointLevel_[pointi] < wantedLevel)
753 dumpCell(mesh_.faceOwner()[facei]);
754 if (mesh_.isInternalFace(facei))
756 dumpCell(mesh_.faceNeighbour()[facei]);
762 <<
" startFp:" << startFp
763 <<
" wantedLevel:" << wantedLevel
766 else if (pointLevel_[pointi] == wantedLevel)
781 dumpCell(mesh_.faceOwner()[facei]);
782 if (mesh_.isInternalFace(facei))
784 dumpCell(mesh_.faceNeighbour()[facei]);
790 <<
" startFp:" << startFp
791 <<
" wantedLevel:" << wantedLevel
801 const face&
f = mesh_.faces()[facei];
805 return pointLevel_[
f[findMaxLevel(
f)]];
809 label ownLevel = cellLevel_[mesh_.faceOwner()[facei]];
811 if (countAnchors(
f, ownLevel) == 4)
815 else if (countAnchors(
f, ownLevel+1) == 4)
827 void Foam::hexRef8::checkInternalOrientation
840 const vector areaNorm(compactFace.areaNormal(compactPoints));
842 const vector dir(neiPt - ownPt);
844 if ((dir & areaNorm) < 0)
847 <<
"cell:" << celli <<
" old face:" << facei
848 <<
" newFace:" << newFace <<
endl
849 <<
" coords:" << compactPoints
850 <<
" ownPt:" << ownPt
851 <<
" neiPt:" << neiPt
855 const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
857 const scalar
s = (fcToOwn & areaNorm) / (dir & areaNorm);
859 if (s < 0.1 || s > 0.9)
862 <<
"cell:" << celli <<
" old face:" << facei
863 <<
" newFace:" << newFace <<
endl
864 <<
" coords:" << compactPoints
865 <<
" ownPt:" << ownPt
866 <<
" neiPt:" << neiPt
873 void Foam::hexRef8::checkBoundaryOrientation
875 polyTopoChange& meshMod,
879 const point& boundaryPt,
883 face compactFace(
identity(newFace.size()));
884 pointField compactPoints(meshMod.points(), newFace);
886 const vector areaNorm(compactFace.areaNormal(compactPoints));
888 const vector dir(boundaryPt - ownPt);
890 if ((dir & areaNorm) < 0)
893 <<
"cell:" << celli <<
" old face:" << facei
894 <<
" newFace:" << newFace
895 <<
" coords:" << compactPoints
896 <<
" ownPt:" << ownPt
897 <<
" boundaryPt:" << boundaryPt
901 const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
903 const scalar
s = (fcToOwn & dir) /
magSqr(dir);
905 if (s < 0.7 || s > 1.3)
908 <<
"cell:" << celli <<
" old face:" << facei
909 <<
" newFace:" << newFace
910 <<
" coords:" << compactPoints
911 <<
" ownPt:" << ownPt
912 <<
" boundaryPt:" << boundaryPt
922 void Foam::hexRef8::insertEdgeSplit
927 DynamicList<label>& verts
930 if (
p0 < mesh_.nPoints() && p1 < mesh_.nPoints())
934 if (edgeI != -1 && edgeMidPoint[edgeI] != -1)
936 verts.append(edgeMidPoint[edgeI]);
950 Foam::label Foam::hexRef8::storeMidPointInfo
958 const bool faceOrder,
959 const label edgeMidPointi,
960 const label anchorPointi,
961 const label faceMidPointi,
963 Map<edge>& midPointToAnchors,
964 Map<edge>& midPointToFaceMids,
965 polyTopoChange& meshMod
970 bool changed =
false;
971 bool haveTwoAnchors =
false;
973 auto edgeMidFnd = midPointToAnchors.find(edgeMidPointi);
975 if (!edgeMidFnd.found())
977 midPointToAnchors.insert(edgeMidPointi, edge(anchorPointi, -1));
981 edge&
e = edgeMidFnd.val();
983 if (anchorPointi !=
e[0])
992 if (
e[0] != -1 &&
e[1] != -1)
994 haveTwoAnchors =
true;
998 bool haveTwoFaceMids =
false;
1000 auto faceMidFnd = midPointToFaceMids.find(edgeMidPointi);
1002 if (!faceMidFnd.found())
1004 midPointToFaceMids.insert(edgeMidPointi, edge(faceMidPointi, -1));
1008 edge&
e = faceMidFnd.val();
1010 if (faceMidPointi !=
e[0])
1014 e[1] = faceMidPointi;
1019 if (
e[0] != -1 &&
e[1] != -1)
1021 haveTwoFaceMids =
true;
1028 if (changed && haveTwoAnchors && haveTwoFaceMids)
1030 const edge& anchors = midPointToAnchors[edgeMidPointi];
1031 const edge& faceMids = midPointToFaceMids[edgeMidPointi];
1033 label otherFaceMidPointi = faceMids.otherVertex(faceMidPointi);
1040 DynamicList<label> newFaceVerts(4);
1041 if (faceOrder == (mesh_.faceOwner()[facei] == celli))
1043 newFaceVerts.append(faceMidPointi);
1054 newFaceVerts.append(edgeMidPointi);
1064 newFaceVerts.append(otherFaceMidPointi);
1065 newFaceVerts.append(cellMidPoint[celli]);
1069 newFaceVerts.append(otherFaceMidPointi);
1079 newFaceVerts.append(edgeMidPointi);
1089 newFaceVerts.append(faceMidPointi);
1090 newFaceVerts.append(cellMidPoint[celli]);
1094 newFace.transfer(newFaceVerts);
1096 label anchorCell0 = getAnchorCell
1104 label anchorCell1 = getAnchorCell
1110 anchors.otherVertex(anchorPointi)
1117 if (anchorCell0 < anchorCell1)
1122 ownPt = mesh_.points()[anchorPointi];
1123 neiPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1132 ownPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1133 neiPt = mesh_.points()[anchorPointi];
1140 if (anchorCell0 < anchorCell1)
1142 ownPt = mesh_.points()[anchorPointi];
1143 neiPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1147 ownPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1148 neiPt = mesh_.points()[anchorPointi];
1151 checkInternalOrientation
1162 return addInternalFace
1180 void Foam::hexRef8::createInternalFaces
1190 polyTopoChange& meshMod
1196 const cell& cFaces = mesh_.cells()[celli];
1197 const label cLevel = cellLevel_[celli];
1200 Map<edge> midPointToAnchors(24);
1202 Map<edge> midPointToFaceMids(24);
1205 DynamicList<label> storage;
1209 label nFacesAdded = 0;
1213 label facei = cFaces[i];
1215 const face&
f = mesh_.faces()[facei];
1216 const labelList& fEdges = mesh_.faceEdges(facei, storage);
1221 label faceMidPointi = -1;
1223 label nAnchors = countAnchors(
f, cLevel);
1231 label anchorFp = -1;
1235 if (pointLevel_[
f[fp]] <= cLevel)
1243 label edgeMid = findLevel
1247 f.fcIndex(anchorFp),
1251 label faceMid = findLevel
1260 faceMidPointi =
f[faceMid];
1262 else if (nAnchors == 4)
1267 faceMidPointi = faceMidPoint[facei];
1271 dumpCell(mesh_.faceOwner()[facei]);
1272 if (mesh_.isInternalFace(facei))
1274 dumpCell(mesh_.faceNeighbour()[facei]);
1278 <<
"nAnchors:" << nAnchors
1279 <<
" facei:" << facei
1291 label point0 =
f[fp0];
1293 if (pointLevel_[point0] <= cLevel)
1302 label edgeMidPointi = -1;
1304 label fp1 =
f.fcIndex(fp0);
1306 if (pointLevel_[
f[fp1]] <= cLevel)
1309 label edgeI = fEdges[fp0];
1311 edgeMidPointi = edgeMidPoint[edgeI];
1313 if (edgeMidPointi == -1)
1317 const labelList& cPoints = mesh_.cellPoints(celli);
1320 <<
"cell:" << celli <<
" cLevel:" << cLevel
1321 <<
" cell points:" << cPoints
1324 <<
" face:" << facei
1328 <<
" faceAnchorLevel:" << faceAnchorLevel[facei]
1329 <<
" faceMidPoint:" << faceMidPoint[facei]
1330 <<
" faceMidPointi:" << faceMidPointi
1338 label edgeMid = findLevel(facei,
f, fp1,
true, cLevel+1);
1340 edgeMidPointi =
f[edgeMid];
1343 label newFacei = storeMidPointInfo
1366 if (nFacesAdded == 12)
1377 label fpMin1 =
f.rcIndex(fp0);
1379 if (pointLevel_[
f[fpMin1]] <= cLevel)
1382 label edgeI = fEdges[fpMin1];
1384 edgeMidPointi = edgeMidPoint[edgeI];
1386 if (edgeMidPointi == -1)
1390 const labelList& cPoints = mesh_.cellPoints(celli);
1393 <<
"cell:" << celli <<
" cLevel:" << cLevel
1394 <<
" cell points:" << cPoints
1397 <<
" face:" << facei
1401 <<
" faceAnchorLevel:" << faceAnchorLevel[facei]
1402 <<
" faceMidPoint:" << faceMidPoint[facei]
1403 <<
" faceMidPointi:" << faceMidPointi
1411 label edgeMid = findLevel
1420 edgeMidPointi =
f[edgeMid];
1423 newFacei = storeMidPointInfo
1446 if (nFacesAdded == 12)
1454 if (nFacesAdded == 12)
1462 void Foam::hexRef8::walkFaceToMid
1467 const label startFp,
1468 DynamicList<label>& faceVerts
1471 const face&
f = mesh_.faces()[facei];
1472 const labelList& fEdges = mesh_.faceEdges(facei);
1481 if (edgeMidPoint[fEdges[fp]] >= 0)
1483 faceVerts.
append(edgeMidPoint[fEdges[fp]]);
1488 if (pointLevel_[
f[fp]] <= cLevel)
1494 else if (pointLevel_[
f[fp]] == cLevel+1)
1501 else if (pointLevel_[
f[fp]] == cLevel+2)
1504 faceVerts.append(
f[fp]);
1511 void Foam::hexRef8::walkFaceFromMid
1516 const label startFp,
1517 DynamicList<label>& faceVerts
1520 const face&
f = mesh_.faces()[facei];
1521 const labelList& fEdges = mesh_.faceEdges(facei);
1523 label fp =
f.rcIndex(startFp);
1527 if (pointLevel_[
f[fp]] <= cLevel)
1532 else if (pointLevel_[
f[fp]] == cLevel+1)
1538 else if (pointLevel_[
f[fp]] == cLevel+2)
1548 if (edgeMidPoint[fEdges[fp]] >= 0)
1550 faceVerts.
append(edgeMidPoint[fEdges[fp]]);
1566 Foam::label Foam::hexRef8::faceConsistentRefinement
1576 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1578 label own = mesh_.faceOwner()[facei];
1579 label ownLevel = cellLevel[own] + refineCell.get(own);
1581 label nei = mesh_.faceNeighbour()[facei];
1582 label neiLevel = cellLevel[nei] + refineCell.get(nei);
1584 if (ownLevel > (neiLevel+1))
1588 refineCell.set(nei);
1592 refineCell.unset(own);
1596 else if (neiLevel > (ownLevel+1))
1600 refineCell.set(own);
1604 refineCell.unset(nei);
1613 labelList neiLevel(mesh_.nBoundaryFaces());
1617 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1619 neiLevel[i] = cellLevel[own] + refineCell.get(own);
1628 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1629 label ownLevel = cellLevel[own] + refineCell.get(own);
1631 if (ownLevel > (neiLevel[i]+1))
1635 refineCell.unset(own);
1639 else if (neiLevel[i] > (ownLevel+1))
1643 refineCell.set(own);
1654 void Foam::hexRef8::checkWantedRefinementLevels
1660 bitSet refineCell(mesh_.nCells(), cellsToRefine);
1662 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1664 label own = mesh_.faceOwner()[facei];
1665 label ownLevel = cellLevel[own] + refineCell.get(own);
1667 label nei = mesh_.faceNeighbour()[facei];
1668 label neiLevel = cellLevel[nei] + refineCell.get(nei);
1670 if (
mag(ownLevel-neiLevel) > 1)
1676 <<
" current level:" << cellLevel[own]
1677 <<
" level after refinement:" << ownLevel
1679 <<
"neighbour cell:" << nei
1680 <<
" current level:" << cellLevel[nei]
1681 <<
" level after refinement:" << neiLevel
1683 <<
"which does not satisfy 2:1 constraints anymore."
1690 labelList neiLevel(mesh_.nBoundaryFaces());
1694 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1696 neiLevel[i] = cellLevel[own] + refineCell.get(own);
1705 label facei = i + mesh_.nInternalFaces();
1707 label own = mesh_.faceOwner()[facei];
1708 label ownLevel = cellLevel[own] + refineCell.get(own);
1710 if (
mag(ownLevel - neiLevel[i]) > 1)
1712 label patchi = mesh_.boundaryMesh().whichPatch(facei);
1716 <<
"Celllevel does not satisfy 2:1 constraint."
1717 <<
" On coupled face "
1719 <<
" on patch " << patchi <<
" "
1720 << mesh_.boundaryMesh()[patchi].name()
1721 <<
" owner cell " << own
1722 <<
" current level:" << cellLevel[own]
1723 <<
" level after refinement:" << ownLevel
1725 <<
" (coupled) neighbour cell will get refinement "
1738 Pout<<
"hexRef8::setInstance(const fileName& inst) : "
1739 <<
"Resetting file instance to " << inst <<
endl;
1742 cellLevel_.instance() = inst;
1743 pointLevel_.instance() = inst;
1744 level0Edge_.instance() = inst;
1745 history_.instance() = inst;
1749 void Foam::hexRef8::collectLevelPoints
1758 if (pointLevel_[
f[fp]] <= level)
1766 void Foam::hexRef8::collectLevelPoints
1771 DynamicList<label>&
points
1776 label pointi = meshPoints[
f[fp]];
1777 if (pointLevel_[pointi] <= level)
1786 bool Foam::hexRef8::matchHexShape
1789 const label cellLevel,
1790 DynamicList<face>& quads
1793 const cell& cFaces = mesh_.cells()[celli];
1796 DynamicList<label> verts(4);
1804 label facei = cFaces[i];
1805 const face&
f = mesh_.faces()[facei];
1808 collectLevelPoints(
f, cellLevel, verts);
1809 if (verts.size() == 4)
1811 if (mesh_.faceOwner()[facei] != celli)
1815 quads.append(face(0));
1822 if (quads.size() < 6)
1824 Map<labelList> pointFaces(2*cFaces.size());
1828 label facei = cFaces[i];
1829 const face&
f = mesh_.faces()[facei];
1837 collectLevelPoints(
f, cellLevel, verts);
1838 if (verts.size() == 1)
1844 label pointi =
f[fp];
1845 if (pointLevel_[pointi] == cellLevel+1)
1847 auto iter = pointFaces.find(pointi);
1852 pFaces.appendUniq(facei);
1878 label facei =
pFaces[pFacei];
1879 const face&
f = mesh_.faces()[facei];
1880 if (mesh_.faceOwner()[facei] == celli)
1882 fourFaces[pFacei] =
f;
1886 fourFaces[pFacei] =
f.reverseFace();
1892 SubList<face>(fourFaces, fourFaces.size()),
1897 if (edgeLoops.size() == 1)
1903 bigFace.meshPoints(),
1904 bigFace.edgeLoops()[0],
1909 if (verts.size() == 4)
1911 quads.append(face(0));
1920 return (quads.size() == 6);
1935 mesh_.facesInstance(),
1948 mesh_.facesInstance(),
1961 mesh_.facesInstance(),
1974 "refinementHistory",
1975 mesh_.facesInstance(),
1982 (readHistory ? mesh_.nCells() : 0)
1984 faceRemover_(mesh_, GREAT),
1985 savedPointLevel_(0),
2000 <<
"History enabled but number of visible cells "
2003 <<
" is not equal to the number of cells in the mesh "
2010 cellLevel_.size() != mesh_.nCells()
2011 || pointLevel_.size() != mesh_.nPoints()
2015 <<
"Restarted from inconsistent cellLevel or pointLevel files."
2019 <<
"Number of cells in mesh:" << mesh_.nCells()
2020 <<
" does not equal size of cellLevel:" << cellLevel_.size() <<
endl
2021 <<
"Number of points in mesh:" << mesh_.nPoints()
2022 <<
" does not equal size of pointLevel:" << pointLevel_.size()
2040 Foam::hexRef8::hexRef8
2046 const scalar level0Edge
2055 mesh_.facesInstance(),
2068 mesh_.facesInstance(),
2081 mesh_.facesInstance(),
2092 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2099 "refinementHistory",
2100 mesh_.facesInstance(),
2108 faceRemover_(mesh_, GREAT),
2109 savedPointLevel_(0),
2112 if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
2115 <<
"History enabled but number of visible cells in it "
2116 << history_.visibleCells().size()
2117 <<
" is not equal to the number of cells in the mesh "
2123 cellLevel_.size() != mesh_.nCells()
2124 || pointLevel_.size() != mesh_.nPoints()
2128 <<
"Incorrect cellLevel or pointLevel size." <<
endl
2129 <<
"Number of cells in mesh:" << mesh_.nCells()
2130 <<
" does not equal size of cellLevel:" << cellLevel_.size() <<
endl
2131 <<
"Number of points in mesh:" << mesh_.nPoints()
2132 <<
" does not equal size of pointLevel:" << pointLevel_.size()
2137 checkRefinementLevels(-1,
labelList(0));
2149 Foam::hexRef8::hexRef8
2154 const scalar level0Edge
2163 mesh_.facesInstance(),
2176 mesh_.facesInstance(),
2189 mesh_.facesInstance(),
2200 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2207 "refinementHistory",
2208 mesh_.facesInstance(),
2218 faceRemover_(mesh_, GREAT),
2219 savedPointLevel_(0),
2224 cellLevel_.size() != mesh_.nCells()
2225 || pointLevel_.size() != mesh_.nPoints()
2229 <<
"Incorrect cellLevel or pointLevel size." <<
endl
2230 <<
"Number of cells in mesh:" << mesh_.nCells()
2231 <<
" does not equal size of cellLevel:" << cellLevel_.size() <<
endl
2232 <<
"Number of points in mesh:" << mesh_.nPoints()
2233 <<
" does not equal size of pointLevel:" << pointLevel_.size()
2238 checkRefinementLevels(-1,
labelList(0));
2267 label nChanged = faceConsistentRefinement
2278 Pout<<
"hexRef8::consistentRefinement : Changed " << nChanged
2279 <<
" refinement levels due to 2:1 conflicts."
2294 checkWantedRefinementLevels(cellLevel, newCellsToRefine);
2297 return newCellsToRefine;
2309 const label maxFaceDiff,
2312 const label maxPointDiff,
2316 const labelList& faceOwner = mesh_.faceOwner();
2317 const labelList& faceNeighbour = mesh_.faceNeighbour();
2320 if (maxFaceDiff <= 0)
2323 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl
2342 forAll(allCellInfo, celli)
2348 maxFaceDiff*(cellLevel_[celli]+1),
2349 maxFaceDiff*cellLevel_[celli]
2356 label celli = cellsToRefine[i];
2358 allCellInfo[celli].count() = allCellInfo[celli].refinementCount();
2368 int dummyTrackData = 0;
2376 label facei = facesToCheck[i];
2378 if (allFaceInfo[facei].valid(dummyTrackData))
2382 <<
"Argument facesToCheck seems to have duplicate entries!"
2384 <<
"face:" << facei <<
" occurs at positions "
2392 if (mesh_.isInternalFace(facei))
2397 const refinementData& neiData = allCellInfo[faceNeighbour[facei]];
2400 label faceRefineCount;
2403 faceCount = neiData.
count() + maxFaceDiff;
2404 faceRefineCount = faceCount + maxFaceDiff;
2408 faceCount = ownData.
count() + maxFaceDiff;
2409 faceRefineCount = faceCount + maxFaceDiff;
2412 seedFaces.append(facei);
2413 seedFacesInfo.append
2421 allFaceInfo[facei] = seedFacesInfo.last();
2425 label faceCount = ownData.
count() + maxFaceDiff;
2426 label faceRefineCount = faceCount + maxFaceDiff;
2428 seedFaces.append(facei);
2429 seedFacesInfo.append
2437 allFaceInfo[facei] = seedFacesInfo.last();
2445 forAll(faceNeighbour, facei)
2448 if (!allFaceInfo[facei].valid(dummyTrackData))
2450 label own = faceOwner[facei];
2451 label nei = faceNeighbour[facei];
2455 if (allCellInfo[own].
count() > allCellInfo[nei].
count())
2457 allFaceInfo[facei].updateFace
2467 seedFacesInfo.append(allFaceInfo[facei]);
2469 else if (allCellInfo[own].
count() < allCellInfo[nei].
count())
2471 allFaceInfo[facei].updateFace
2480 seedFaces.append(facei);
2481 seedFacesInfo.append(allFaceInfo[facei]);
2489 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
2492 if (!allFaceInfo[facei].valid(dummyTrackData))
2494 label own = faceOwner[facei];
2507 seedFaces.append(facei);
2508 seedFacesInfo.append(faceData);
2526 Pout<<
"hexRef8::consistentSlowRefinement : Seeded "
2527 << seedFaces.size() <<
" faces between cells with different"
2528 <<
" refinement level." <<
endl;
2532 levelCalc.
setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
2534 seedFacesInfo.clear();
2537 levelCalc.
iterate(mesh_.globalData().nTotalFaces()+1);
2545 if (maxPointDiff == -1)
2555 forAll(maxPointCount, pointi)
2557 label& pLevel = maxPointCount[pointi];
2559 const labelList& pCells = mesh_.pointCells(pointi);
2563 pLevel =
max(pLevel, allCellInfo[pCells[i]].
count());
2584 label pointi = pointsToCheck[i];
2589 const labelList& pCells = mesh_.pointCells(pointi);
2593 label celli = pCells[pCelli];
2601 maxPointCount[pointi]
2602 >
cellInfo.count() + maxFaceDiff*maxPointDiff
2610 const cell& cFaces = mesh_.cells()[celli];
2614 label facei = cFaces[cFacei];
2627 if (faceData.
count() > allFaceInfo[facei].count())
2629 changedFacesInfo.insert(facei, faceData);
2636 label nChanged = changedFacesInfo.size();
2646 seedFaces.setCapacity(changedFacesInfo.size());
2647 seedFacesInfo.setCapacity(changedFacesInfo.size());
2651 seedFaces.append(iter.key());
2652 seedFacesInfo.append(iter.val());
2659 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2661 label own = mesh_.faceOwner()[facei];
2664 + (allCellInfo[own].isRefined() ? 1 : 0);
2666 label nei = mesh_.faceNeighbour()[facei];
2669 + (allCellInfo[nei].isRefined() ? 1 : 0);
2671 if (
mag(ownLevel-neiLevel) > 1)
2677 <<
" current level:" << cellLevel_[own]
2678 <<
" current refData:" << allCellInfo[own]
2679 <<
" level after refinement:" << ownLevel
2681 <<
"neighbour cell:" << nei
2682 <<
" current level:" << cellLevel_[nei]
2683 <<
" current refData:" << allCellInfo[nei]
2684 <<
" level after refinement:" << neiLevel
2686 <<
"which does not satisfy 2:1 constraints anymore." <<
nl
2687 <<
"face:" << facei <<
" faceRefData:" << allFaceInfo[facei]
2696 labelList neiLevel(mesh_.nBoundaryFaces());
2697 labelList neiCount(mesh_.nBoundaryFaces());
2698 labelList neiRefCount(mesh_.nBoundaryFaces());
2702 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
2703 neiLevel[i] = cellLevel_[own];
2704 neiCount[i] = allCellInfo[own].count();
2705 neiRefCount[i] = allCellInfo[own].refinementCount();
2716 label facei = i+mesh_.nInternalFaces();
2718 label own = mesh_.faceOwner()[facei];
2721 + (allCellInfo[own].isRefined() ? 1 : 0);
2725 + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
2727 if (
mag(ownLevel - nbrLevel) > 1)
2730 label patchi = mesh_.boundaryMesh().whichPatch(facei);
2733 <<
"Celllevel does not satisfy 2:1 constraint."
2734 <<
" On coupled face "
2736 <<
" refData:" << allFaceInfo[facei]
2737 <<
" on patch " << patchi <<
" "
2738 << mesh_.boundaryMesh()[patchi].name() <<
nl
2739 <<
"owner cell " << own
2740 <<
" current level:" << cellLevel_[own]
2741 <<
" current count:" << allCellInfo[own].count()
2742 <<
" current refCount:"
2743 << allCellInfo[own].refinementCount()
2744 <<
" level after refinement:" << ownLevel
2746 <<
"(coupled) neighbour cell"
2747 <<
" has current level:" << neiLevel[i]
2748 <<
" current count:" << neiCount[i]
2749 <<
" current refCount:" << neiRefCount[i]
2750 <<
" level after refinement:" << nbrLevel
2760 forAll(allCellInfo, celli)
2762 if (allCellInfo[celli].isRefined())
2772 forAll(allCellInfo, celli)
2774 if (allCellInfo[celli].isRefined())
2776 newCellsToRefine[nRefined++] = celli;
2782 Pout<<
"hexRef8::consistentSlowRefinement : From "
2783 << cellsToRefine.size() <<
" to " << newCellsToRefine.size()
2784 <<
" cells to refine." <<
endl;
2787 return newCellsToRefine;
2793 const label maxFaceDiff,
2798 const labelList& faceOwner = mesh_.faceOwner();
2799 const labelList& faceNeighbour = mesh_.faceNeighbour();
2801 if (maxFaceDiff <= 0)
2804 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl
2808 const scalar level0Size = 2*maxFaceDiff*level0EdgeLength();
2823 int dummyTrackData = 0;
2829 label celli = cellsToRefine[i];
2834 mesh_.cellCentres()[celli],
2839 forAll(allCellInfo, celli)
2841 if (!allCellInfo[celli].valid(dummyTrackData))
2846 mesh_.cellCentres()[celli],
2862 label facei = facesToCheck[i];
2864 if (allFaceInfo[facei].valid(dummyTrackData))
2868 <<
"Argument facesToCheck seems to have duplicate entries!"
2870 <<
"face:" << facei <<
" occurs at positions "
2875 label own = faceOwner[facei];
2879 allCellInfo[own].valid(dummyTrackData)
2880 ? allCellInfo[own].originLevel()
2884 if (!mesh_.isInternalFace(facei))
2888 const point& fc = mesh_.faceCentres()[facei];
2897 allFaceInfo[facei].updateFace
2909 label nei = faceNeighbour[facei];
2913 allCellInfo[nei].valid(dummyTrackData)
2914 ? allCellInfo[nei].originLevel()
2918 if (ownLevel == neiLevel)
2921 allFaceInfo[facei].updateFace
2930 allFaceInfo[facei].updateFace
2943 allFaceInfo[facei].updateFace
2952 allFaceInfo[facei].updateFace
2964 seedFacesInfo.append(allFaceInfo[facei]);
2971 forAll(faceNeighbour, facei)
2974 if (!allFaceInfo[facei].valid(dummyTrackData))
2976 label own = faceOwner[facei];
2980 allCellInfo[own].valid(dummyTrackData)
2981 ? allCellInfo[own].originLevel()
2985 label nei = faceNeighbour[facei];
2989 allCellInfo[nei].valid(dummyTrackData)
2990 ? allCellInfo[nei].originLevel()
2994 if (ownLevel > neiLevel)
2998 allFaceInfo[facei].updateFace
3007 seedFacesInfo.append(allFaceInfo[facei]);
3009 else if (neiLevel > ownLevel)
3011 seedFaces.append(facei);
3012 allFaceInfo[facei].updateFace
3021 seedFacesInfo.append(allFaceInfo[facei]);
3027 seedFacesInfo.shrink();
3037 mesh_.globalData().nTotalCells()+1,
3078 label celli = cellsToRefine[i];
3080 allCellInfo[celli].originLevel() =
sizeof(label)*8-2;
3081 allCellInfo[celli].origin() = cc[celli];
3088 forAll(allCellInfo, celli)
3090 label wanted = allCellInfo[celli].wantedLevel(cc[celli]);
3092 if (wanted > cellLevel_[celli]+1)
3097 faceConsistentRefinement(
true, cellLevel_,
refineCell);
3101 label nChanged = faceConsistentRefinement(
true, cellLevel_,
refineCell);
3107 Pout<<
"hexRef8::consistentSlowRefinement2 : Changed " << nChanged
3108 <<
" refinement levels due to 2:1 conflicts."
3123 Pout<<
"hexRef8::consistentSlowRefinement2 : From "
3124 << cellsToRefine.size() <<
" to " << newCellsToRefine.size()
3125 <<
" cells to refine." <<
endl;
3130 cellSet cellsIn(mesh_,
"cellsToRefineIn", cellsToRefine);
3131 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
3132 << cellsIn.size() <<
" to cellSet "
3137 cellSet cellsOut(mesh_,
"cellsToRefineOut", newCellsToRefine);
3138 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
3139 << cellsOut.size() <<
" to cellSet "
3148 label nChanged = faceConsistentRefinement(
true, cellLevel_,
refineCell);
3153 mesh_,
"cellsToRefineOut2", newCellsToRefine.size()
3159 cellsOut2.insert(celli);
3162 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
3163 << cellsOut2.size() <<
" to cellSet "
3164 << cellsOut2.objectPath() <<
endl;
3172 if (
refineCell.test(celli) && !savedRefineCell.test(celli))
3176 <<
"Cell:" << celli <<
" cc:"
3177 << mesh_.cellCentres()[celli]
3178 <<
" was not marked for refinement but does not obey"
3179 <<
" 2:1 constraints."
3186 return newCellsToRefine;
3199 Pout<<
"hexRef8::setRefinement :"
3200 <<
" Checking initial mesh just to make sure" <<
endl;
3209 savedPointLevel_.clear();
3210 savedCellLevel_.clear();
3215 forAll(cellLevel_, celli)
3217 newCellLevel.append(cellLevel_[celli]);
3220 forAll(pointLevel_, pointi)
3222 newPointLevel.append(pointLevel_[pointi]);
3228 Pout<<
"hexRef8::setRefinement :"
3229 <<
" Allocating " << cellLabels.size() <<
" cell midpoints."
3237 labelList cellMidPoint(mesh_.nCells(), -1);
3241 label celli = cellLabels[i];
3243 label anchorPointi = mesh_.faces()[mesh_.cells()[celli][0]][0];
3249 mesh_.cellCentres()[celli],
3256 newPointLevel(cellMidPoint[celli]) = cellLevel_[celli]+1;
3262 cellSet splitCells(mesh_,
"splitCells", cellLabels.size());
3264 forAll(cellMidPoint, celli)
3266 if (cellMidPoint[celli] >= 0)
3268 splitCells.insert(celli);
3272 Pout<<
"hexRef8::setRefinement : Dumping " << splitCells.size()
3273 <<
" cells to split to cellSet " << splitCells.objectPath()
3286 Pout<<
"hexRef8::setRefinement :"
3287 <<
" Allocating edge midpoints."
3296 labelList edgeMidPoint(mesh_.nEdges(), -1);
3299 forAll(cellMidPoint, celli)
3301 if (cellMidPoint[celli] >= 0)
3303 const labelList& cEdges = mesh_.cellEdges(celli);
3307 label edgeI = cEdges[i];
3309 const edge&
e = mesh_.edges()[edgeI];
3313 pointLevel_[
e[0]] <= cellLevel_[celli]
3314 && pointLevel_[
e[1]] <= cellLevel_[celli]
3317 edgeMidPoint[edgeI] = 12345;
3344 forAll(edgeMidPoint, edgeI)
3346 if (edgeMidPoint[edgeI] >= 0)
3349 edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points());
3357 point(-GREAT, -GREAT, -GREAT)
3362 forAll(edgeMidPoint, edgeI)
3364 if (edgeMidPoint[edgeI] >= 0)
3369 const edge&
e = mesh_.edges()[edgeI];
3382 newPointLevel(edgeMidPoint[edgeI]) =
3395 OFstream str(mesh_.time().path()/
"edgeMidPoint.obj");
3397 forAll(edgeMidPoint, edgeI)
3399 if (edgeMidPoint[edgeI] >= 0)
3401 const edge&
e = mesh_.edges()[edgeI];
3407 Pout<<
"hexRef8::setRefinement :"
3408 <<
" Dumping edge centres to split to file " << str.
name() <<
endl;
3418 Pout<<
"hexRef8::setRefinement :"
3419 <<
" Allocating face midpoints."
3425 labelList faceAnchorLevel(mesh_.nFaces());
3427 for (label facei = 0; facei < mesh_.nFaces(); facei++)
3429 faceAnchorLevel[facei] = faceLevel(facei);
3434 labelList faceMidPoint(mesh_.nFaces(), -1);
3439 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
3441 if (faceAnchorLevel[facei] >= 0)
3443 label own = mesh_.faceOwner()[facei];
3444 label ownLevel = cellLevel_[own];
3445 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3447 label nei = mesh_.faceNeighbour()[facei];
3448 label neiLevel = cellLevel_[nei];
3449 label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
3453 newOwnLevel > faceAnchorLevel[facei]
3454 || newNeiLevel > faceAnchorLevel[facei]
3457 faceMidPoint[facei] = 12345;
3470 labelList newNeiLevel(mesh_.nBoundaryFaces());
3474 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
3475 label ownLevel = cellLevel_[own];
3476 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3478 newNeiLevel[i] = newOwnLevel;
3488 label facei = i+mesh_.nInternalFaces();
3490 if (faceAnchorLevel[facei] >= 0)
3492 label own = mesh_.faceOwner()[facei];
3493 label ownLevel = cellLevel_[own];
3494 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3498 newOwnLevel > faceAnchorLevel[facei]
3499 || newNeiLevel[i] > faceAnchorLevel[facei]
3502 faceMidPoint[facei] = 12345;
3527 mesh_.nBoundaryFaces(),
3528 point(-GREAT, -GREAT, -GREAT)
3533 label facei = i+mesh_.nInternalFaces();
3535 if (faceMidPoint[facei] >= 0)
3537 bFaceMids[i] = mesh_.faceCentres()[facei];
3547 forAll(faceMidPoint, facei)
3549 if (faceMidPoint[facei] >= 0)
3554 const face&
f = mesh_.faces()[facei];
3561 facei < mesh_.nInternalFaces()
3562 ? mesh_.faceCentres()[facei]
3563 : bFaceMids[facei-mesh_.nInternalFaces()]
3573 newPointLevel(faceMidPoint[facei]) = faceAnchorLevel[facei]+1;
3580 faceSet splitFaces(mesh_,
"splitFaces", cellLabels.size());
3582 forAll(faceMidPoint, facei)
3584 if (faceMidPoint[facei] >= 0)
3586 splitFaces.insert(facei);
3590 Pout<<
"hexRef8::setRefinement : Dumping " << splitFaces.size()
3591 <<
" faces to split to faceSet " << splitFaces.objectPath() <<
endl;
3612 Pout<<
"hexRef8::setRefinement :"
3613 <<
" Finding cell anchorPoints (8 per cell)"
3626 forAll(cellMidPoint, celli)
3628 if (cellMidPoint[celli] >= 0)
3630 cellAnchorPoints[celli].setSize(8);
3634 forAll(pointLevel_, pointi)
3636 const labelList& pCells = mesh_.pointCells(pointi);
3640 label celli = pCells[pCelli];
3644 cellMidPoint[celli] >= 0
3645 && pointLevel_[pointi] <= cellLevel_[celli]
3648 if (nAnchorPoints[celli] == 8)
3653 <<
" of level " << cellLevel_[celli]
3654 <<
" uses more than 8 points of equal or"
3655 <<
" lower level" <<
nl
3656 <<
"Points so far:" << cellAnchorPoints[celli]
3660 cellAnchorPoints[celli][nAnchorPoints[celli]++] = pointi;
3665 forAll(cellMidPoint, celli)
3667 if (cellMidPoint[celli] >= 0)
3669 if (nAnchorPoints[celli] != 8)
3673 const labelList& cPoints = mesh_.cellPoints(celli);
3677 <<
" of level " << cellLevel_[celli]
3678 <<
" does not seem to have 8 points of equal or"
3679 <<
" lower level" <<
endl
3680 <<
"cellPoints:" << cPoints <<
endl
3695 Pout<<
"hexRef8::setRefinement :"
3696 <<
" Adding cells (1 per anchorPoint)"
3703 forAll(cellAnchorPoints, celli)
3705 const labelList& cAnchors = cellAnchorPoints[celli];
3707 if (cAnchors.size() == 8)
3709 labelList& cAdded = cellAddedCells[celli];
3716 newCellLevel[celli] = cellLevel_[celli]+1;
3719 for (label i = 1; i < 8; i++)
3729 mesh_.cellZones().whichZone(celli)
3733 newCellLevel(cAdded[i]) = cellLevel_[celli]+1;
3748 Pout<<
"hexRef8::setRefinement :"
3749 <<
" Marking faces to be handled"
3754 bitSet affectedFace(mesh_.nFaces());
3757 forAll(cellMidPoint, celli)
3759 if (cellMidPoint[celli] >= 0)
3761 const cell& cFaces = mesh_.cells()[celli];
3763 affectedFace.
set(cFaces);
3767 forAll(faceMidPoint, facei)
3769 if (faceMidPoint[facei] >= 0)
3771 affectedFace.set(facei);
3775 forAll(edgeMidPoint, edgeI)
3777 if (edgeMidPoint[edgeI] >= 0)
3779 const labelList& eFaces = mesh_.edgeFaces(edgeI);
3781 affectedFace.
set(eFaces);
3792 Pout<<
"hexRef8::setRefinement : Splitting faces" <<
endl;
3795 forAll(faceMidPoint, facei)
3797 if (faceMidPoint[facei] >= 0 && affectedFace.test(facei))
3803 const face&
f = mesh_.faces()[facei];
3807 bool modifiedFace =
false;
3809 label anchorLevel = faceAnchorLevel[facei];
3815 label pointi =
f[fp];
3817 if (pointLevel_[pointi] <= anchorLevel)
3823 faceVerts.
append(pointi);
3839 faceVerts.
append(faceMidPoint[facei]);
3872 if (mesh_.isInternalFace(facei))
3874 label oldOwn = mesh_.faceOwner()[facei];
3875 label oldNei = mesh_.faceNeighbour()[facei];
3877 checkInternalOrientation
3882 mesh_.cellCentres()[oldOwn],
3883 mesh_.cellCentres()[oldNei],
3889 label oldOwn = mesh_.faceOwner()[facei];
3891 checkBoundaryOrientation
3896 mesh_.cellCentres()[oldOwn],
3897 mesh_.faceCentres()[facei],
3906 modifiedFace =
true;
3908 modFace(meshMod, facei, newFace, own, nei);
3912 addFace(meshMod, facei, newFace, own, nei);
3918 affectedFace.unset(facei);
3928 Pout<<
"hexRef8::setRefinement :"
3929 <<
" Adding edge splits to unsplit faces"
3936 forAll(edgeMidPoint, edgeI)
3938 if (edgeMidPoint[edgeI] >= 0)
3942 const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
3946 label facei = eFaces[i];
3948 if (faceMidPoint[facei] < 0 && affectedFace.test(facei))
3952 const face&
f = mesh_.faces()[facei];
3953 const labelList& fEdges = mesh_.faceEdges
3963 newFaceVerts.append(
f[fp]);
3965 label edgeI = fEdges[fp];
3967 if (edgeMidPoint[edgeI] >= 0)
3969 newFaceVerts.
append(edgeMidPoint[edgeI]);
3978 label anchorFp = findMinLevel(
f);
3995 if (mesh_.isInternalFace(facei))
3997 label oldOwn = mesh_.faceOwner()[facei];
3998 label oldNei = mesh_.faceNeighbour()[facei];
4000 checkInternalOrientation
4005 mesh_.cellCentres()[oldOwn],
4006 mesh_.cellCentres()[oldNei],
4012 label oldOwn = mesh_.faceOwner()[facei];
4014 checkBoundaryOrientation
4019 mesh_.cellCentres()[oldOwn],
4020 mesh_.faceCentres()[facei],
4026 modFace(meshMod, facei, newFace, own, nei);
4029 affectedFace.unset(facei);
4041 Pout<<
"hexRef8::setRefinement :"
4042 <<
" Changing owner/neighbour for otherwise unaffected faces"
4046 forAll(affectedFace, facei)
4048 if (affectedFace.test(facei))
4050 const face&
f = mesh_.faces()[facei];
4054 label anchorFp = findMinLevel(
f);
4068 modFace(meshMod, facei,
f, own, nei);
4071 affectedFace.unset(facei);
4086 Pout<<
"hexRef8::setRefinement :"
4087 <<
" Create new internal faces for split cells"
4091 forAll(cellMidPoint, celli)
4093 if (cellMidPoint[celli] >= 0)
4118 forAll(cellMidPoint, celli)
4120 if (cellMidPoint[celli] >= 0)
4122 minPointi =
min(minPointi, cellMidPoint[celli]);
4123 maxPointi =
max(maxPointi, cellMidPoint[celli]);
4126 forAll(faceMidPoint, facei)
4128 if (faceMidPoint[facei] >= 0)
4130 minPointi =
min(minPointi, faceMidPoint[facei]);
4131 maxPointi =
max(maxPointi, faceMidPoint[facei]);
4134 forAll(edgeMidPoint, edgeI)
4136 if (edgeMidPoint[edgeI] >= 0)
4138 minPointi =
min(minPointi, edgeMidPoint[edgeI]);
4139 maxPointi =
max(maxPointi, edgeMidPoint[edgeI]);
4143 if (minPointi !=
labelMax && minPointi != mesh_.nPoints())
4146 <<
"Added point labels not consecutive to existing mesh points."
4148 <<
"mesh_.nPoints():" << mesh_.nPoints()
4149 <<
" minPointi:" << minPointi
4150 <<
" maxPointi:" << maxPointi
4155 pointLevel_.transfer(newPointLevel);
4156 cellLevel_.transfer(newCellLevel);
4159 setInstance(mesh_.facesInstance());
4166 if (history_.active())
4170 Pout<<
"hexRef8::setRefinement :"
4171 <<
" Updating refinement history to " << cellLevel_.size()
4172 <<
" cells" <<
endl;
4176 history_.resize(cellLevel_.size());
4178 forAll(cellAddedCells, celli)
4180 const labelList& addedCells = cellAddedCells[celli];
4182 if (addedCells.size())
4185 history_.storeSplit(celli, addedCells);
4196 label celli = cellLabels[i];
4198 refinedCells[i].
transfer(cellAddedCells[celli]);
4201 return refinedCells;
4212 savedPointLevel_.clear();
4213 savedPointLevel_.resize(2*pointsToStore.size());
4216 label pointi = pointsToStore[i];
4217 savedPointLevel_.insert(pointi, pointLevel_[pointi]);
4220 savedCellLevel_.
clear();
4221 savedCellLevel_.resize(2*cellsToStore.size());
4224 label celli = cellsToStore[i];
4225 savedCellLevel_.insert(celli, cellLevel_[celli]);
4237 updateMesh(map, dummyMap, dummyMap, dummyMap);
4255 Pout<<
"hexRef8::updateMesh :"
4256 <<
" Updating various lists"
4265 Pout<<
"hexRef8::updateMesh :"
4267 <<
" cellMap:" << map.
cellMap().size()
4268 <<
" nCells:" << mesh_.nCells()
4270 <<
" cellLevel_:" << cellLevel_.size()
4272 <<
" pointMap:" << map.
pointMap().size()
4273 <<
" nPoints:" << mesh_.nPoints()
4275 <<
" pointLevel_:" << pointLevel_.size()
4279 if (reverseCellMap.size() == cellLevel_.size())
4285 reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
4293 forAll(cellMap, newCelli)
4295 label oldCelli = cellMap[newCelli];
4299 newCellLevel[newCelli] = -1;
4303 newCellLevel[newCelli] = cellLevel_[oldCelli];
4313 const label newCelli = iter.key();
4314 const label storedCelli = iter.val();
4316 const auto fnd = savedCellLevel_.cfind(storedCelli);
4321 <<
"Problem : trying to restore old value for new cell "
4322 << newCelli <<
" but cannot find old cell " << storedCelli
4323 <<
" in map of stored values " << savedCellLevel_
4326 cellLevel_[newCelli] = fnd.val();
4347 if (reversePointMap.size() == pointLevel_.size())
4350 reorder(reversePointMap, mesh_.nPoints(), -1, pointLevel_);
4357 labelList newPointLevel(pointMap.size());
4361 const label oldPointi = pointMap[
newPointi];
4363 if (oldPointi == -1)
4376 newPointLevel[
newPointi] = pointLevel_[oldPointi];
4379 pointLevel_.
transfer(newPointLevel);
4387 const label storedPointi = iter.val();
4389 auto fnd = savedPointLevel_.find(storedPointi);
4394 <<
"Problem : trying to restore old value for new point "
4395 <<
newPointi <<
" but cannot find old point "
4397 <<
" in map of stored values " << savedPointLevel_
4417 if (history_.active())
4419 history_.updateMesh(map);
4423 setInstance(mesh_.facesInstance());
4426 faceRemover_.updateMesh(map);
4429 cellShapesPtr_.clear();
4444 Pout<<
"hexRef8::subset :"
4445 <<
" Updating various lists"
4449 if (history_.active())
4452 <<
"Subsetting will not work in combination with unrefinement."
4454 <<
"Proceed at your own risk." <<
endl;
4462 forAll(cellMap, newCelli)
4464 newCellLevel[newCelli] = cellLevel_[cellMap[newCelli]];
4467 cellLevel_.transfer(newCellLevel);
4469 if (cellLevel_.found(-1))
4473 <<
"cellLevel_ contains illegal value -1 after mapping:"
4481 labelList newPointLevel(pointMap.size());
4488 pointLevel_.transfer(newPointLevel);
4490 if (pointLevel_.found(-1))
4494 <<
"pointLevel_ contains illegal value -1 after mapping:"
4501 if (history_.active())
4503 history_.subset(pointMap,
faceMap, cellMap);
4507 setInstance(mesh_.facesInstance());
4513 cellShapesPtr_.clear();
4522 Pout<<
"hexRef8::distribute :"
4523 <<
" Updating various lists"
4533 if (history_.active())
4535 history_.distribute(map);
4539 faceRemover_.distribute(map);
4542 cellShapesPtr_.clear();
4548 const scalar smallDim = 1
e-6 * mesh_.bounds().mag();
4552 Pout<<
"hexRef8::checkMesh : Using matching tolerance smallDim:"
4553 << smallDim <<
endl;
4566 nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4584 label facei = pp.
start();
4588 label own = mesh_.faceOwner()[facei];
4589 label bFacei = facei-mesh_.nInternalFaces();
4591 if (!cellToFace.insert(
labelPair(own, nei[bFacei]), facei))
4595 <<
"Faces do not seem to be correct across coupled"
4596 <<
" boundaries" <<
endl
4597 <<
"Coupled face " << facei
4598 <<
" between owner " << own
4599 <<
" on patch " << pp.
name()
4600 <<
" and coupled neighbour " << nei[bFacei]
4601 <<
" has two faces connected to it:"
4603 << cellToFace[
labelPair(own, nei[bFacei])]
4620 neiFaceAreas[i] =
mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
4628 label facei = i+mesh_.nInternalFaces();
4630 const scalar magArea =
mag(mesh_.faceAreas()[facei]);
4632 if (
mag(magArea - neiFaceAreas[i]) > smallDim)
4634 const face&
f = mesh_.faces()[facei];
4635 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4637 dumpCell(mesh_.faceOwner()[facei]);
4640 <<
"Faces do not seem to be correct across coupled"
4641 <<
" boundaries" <<
endl
4642 <<
"Coupled face " << facei
4643 <<
" on patch " << patchi
4644 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4646 <<
" has face area:" << magArea
4647 <<
" (coupled) neighbour face area differs:"
4649 <<
" to within tolerance " << smallDim
4659 labelList nVerts(mesh_.nBoundaryFaces());
4663 nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size();
4671 label facei = i+mesh_.nInternalFaces();
4673 const face&
f = mesh_.faces()[facei];
4675 if (
f.size() != nVerts[i])
4677 dumpCell(mesh_.faceOwner()[facei]);
4679 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4682 <<
"Faces do not seem to be correct across coupled"
4683 <<
" boundaries" <<
endl
4684 <<
"Coupled face " << facei
4685 <<
" on patch " << patchi
4686 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4688 <<
" has size:" <<
f.size()
4689 <<
" (coupled) neighbour face has size:"
4701 pointField anchorPoints(mesh_.nBoundaryFaces());
4705 label facei = i+mesh_.nInternalFaces();
4706 const point& fc = mesh_.faceCentres()[facei];
4707 const face&
f = mesh_.faces()[facei];
4708 const vector anchorVec(mesh_.points()[
f[0]] - fc);
4710 anchorPoints[i] = anchorVec;
4720 label facei = i+mesh_.nInternalFaces();
4721 const point& fc = mesh_.faceCentres()[facei];
4722 const face&
f = mesh_.faces()[facei];
4723 const vector anchorVec(mesh_.points()[
f[0]] - fc);
4725 if (
mag(anchorVec - anchorPoints[i]) > smallDim)
4727 dumpCell(mesh_.faceOwner()[facei]);
4729 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4732 <<
"Faces do not seem to be correct across coupled"
4733 <<
" boundaries" <<
endl
4734 <<
"Coupled face " << facei
4735 <<
" on patch " << patchi
4736 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4738 <<
" has anchor vector:" << anchorVec
4739 <<
" (coupled) neighbour face anchor vector differs:"
4741 <<
" to within tolerance " << smallDim
4749 Pout<<
"hexRef8::checkMesh : Returning" <<
endl;
4756 const label maxPointDiff,
4762 Pout<<
"hexRef8::checkRefinementLevels :"
4763 <<
" Checking 2:1 refinement level" <<
endl;
4768 cellLevel_.size() != mesh_.nCells()
4769 || pointLevel_.size() != mesh_.nPoints()
4773 <<
"cellLevel size should be number of cells"
4774 <<
" and pointLevel size should be number of points."<<
nl
4775 <<
"cellLevel:" << cellLevel_.size()
4776 <<
" mesh.nCells():" << mesh_.nCells() <<
nl
4777 <<
"pointLevel:" << pointLevel_.size()
4778 <<
" mesh.nPoints():" << mesh_.nPoints()
4788 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
4790 label own = mesh_.faceOwner()[facei];
4791 label nei = mesh_.faceNeighbour()[facei];
4793 if (
mag(cellLevel_[own] - cellLevel_[nei]) > 1)
4799 <<
"Celllevel does not satisfy 2:1 constraint." <<
nl
4800 <<
"On face " << facei <<
" owner cell " << own
4801 <<
" has refinement " << cellLevel_[own]
4802 <<
" neighbour cell " << nei <<
" has refinement "
4809 labelList neiLevel(mesh_.nBoundaryFaces());
4813 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4815 neiLevel[i] = cellLevel_[own];
4823 label facei = i+mesh_.nInternalFaces();
4825 label own = mesh_.faceOwner()[facei];
4827 if (
mag(cellLevel_[own] - neiLevel[i]) > 1)
4831 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4834 <<
"Celllevel does not satisfy 2:1 constraint."
4835 <<
" On coupled face " << facei
4836 <<
" on patch " << patchi <<
" "
4837 << mesh_.boundaryMesh()[patchi].name()
4838 <<
" owner cell " << own <<
" has refinement "
4840 <<
" (coupled) neighbour cell has refinement "
4863 forAll(syncPointLevel, pointi)
4865 if (pointLevel_[pointi] != syncPointLevel[pointi])
4868 <<
"PointLevel is not consistent across coupled patches."
4870 <<
"point:" << pointi <<
" coord:" << mesh_.points()[pointi]
4871 <<
" has level " << pointLevel_[pointi]
4872 <<
" whereas the coupled point has level "
4873 << syncPointLevel[pointi]
4881 if (maxPointDiff != -1)
4886 forAll(maxPointLevel, pointi)
4888 const labelList& pCells = mesh_.pointCells(pointi);
4890 label& pLevel = maxPointLevel[pointi];
4894 pLevel =
max(pLevel, cellLevel_[pCells[i]]);
4910 label pointi = pointsToCheck[i];
4912 const labelList& pCells = mesh_.pointCells(pointi);
4916 label celli = pCells[i];
4920 mag(cellLevel_[celli]-maxPointLevel[pointi])
4927 <<
"Too big a difference between"
4928 <<
" point-connected cells." <<
nl
4930 <<
" cellLevel:" << cellLevel_[celli]
4931 <<
" uses point:" << pointi
4932 <<
" coord:" << mesh_.points()[pointi]
4933 <<
" which is also used by a cell with level:"
4934 << maxPointLevel[pointi]
5009 if (!cellShapesPtr_)
5013 Pout<<
"hexRef8::cellShapes() : calculating splitHex cellShapes."
5014 <<
" cellLevel:" << cellLevel_.size()
5015 <<
" pointLevel:" << pointLevel_.size()
5022 label nSplitHex = 0;
5023 label nUnrecognised = 0;
5025 forAll(cellLevel_, celli)
5027 if (meshShapes[celli].model().index() == 0)
5029 label level = cellLevel_[celli];
5033 bool haveQuads = matchHexShape
5043 cellShapesPtr_()[celli] = degenerateMatcher::match(faces);
5054 Pout<<
"hexRef8::cellShapes() :"
5055 <<
" nCells:" << mesh_.nCells() <<
" of which" <<
nl
5056 <<
" primitive:" << (mesh_.nCells()-nSplitHex-nUnrecognised)
5058 <<
" split-hex:" << nSplitHex <<
nl
5059 <<
" poly :" << nUnrecognised <<
nl
5063 return *cellShapesPtr_;
5071 checkRefinementLevels(-1,
labelList(0));
5076 Pout<<
"hexRef8::getSplitPoints :"
5077 <<
" Calculating unrefineable points" <<
endl;
5081 if (!history_.active())
5084 <<
"Only call if constructed with history capability"
5092 labelList splitMaster(mesh_.nPoints(), -1);
5098 for (label pointi = 0; pointi < mesh_.nPoints(); pointi++)
5100 const labelList& pCells = mesh_.pointCells(pointi);
5102 if (pCells.size() != 8)
5104 splitMaster[pointi] = -2;
5109 const labelList& visibleCells = history_.visibleCells();
5111 forAll(visibleCells, celli)
5113 const labelList& cPoints = mesh_.cellPoints(celli);
5115 if (visibleCells[celli] != -1 && history_.parentIndex(celli) >= 0)
5117 label parentIndex = history_.parentIndex(celli);
5122 label pointi = cPoints[i];
5124 label masterCelli = splitMaster[pointi];
5126 if (masterCelli == -1)
5133 splitMaster[pointi] = parentIndex;
5134 splitMasterLevel[pointi] = cellLevel_[celli] - 1;
5136 else if (masterCelli == -2)
5142 (masterCelli != parentIndex)
5143 || (splitMasterLevel[pointi] != cellLevel_[celli] - 1)
5148 splitMaster[pointi] = -2;
5157 label pointi = cPoints[i];
5159 splitMaster[pointi] = -2;
5167 label facei = mesh_.nInternalFaces();
5168 facei < mesh_.nFaces();
5172 const face&
f = mesh_.faces()[facei];
5176 splitMaster[
f[fp]] = -2;
5183 label nSplitPoints = 0;
5185 forAll(splitMaster, pointi)
5187 if (splitMaster[pointi] >= 0)
5196 forAll(splitMaster, pointi)
5198 if (splitMaster[pointi] >= 0)
5200 splitPoints[nSplitPoints++] = pointi;
5278 Pout<<
"hexRef8::consistentUnrefinement :"
5279 <<
" Determining 2:1 consistent unrefinement" <<
endl;
5285 <<
"maxSet not implemented yet."
5295 bitSet unrefinePoint(mesh_.nPoints(), pointsToUnrefine);
5302 bitSet unrefineCell(mesh_.nCells());
5304 forAll(unrefinePoint, pointi)
5306 if (unrefinePoint.test(pointi))
5308 const labelList& pCells = mesh_.pointCells(pointi);
5310 unrefineCell.
set(pCells);
5322 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
5324 label own = mesh_.faceOwner()[facei];
5325 label nei = mesh_.faceNeighbour()[facei];
5327 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5328 label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
5330 if (ownLevel < (neiLevel-1))
5337 unrefineCell.set(nei);
5348 if (!unrefineCell.test(own))
5354 unrefineCell.unset(own);
5358 else if (neiLevel < (ownLevel-1))
5362 unrefineCell.set(own);
5366 if (!unrefineCell.test(nei))
5372 unrefineCell.unset(nei);
5380 labelList neiLevel(mesh_.nBoundaryFaces());
5384 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
5386 neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
5394 label facei = i+mesh_.nInternalFaces();
5395 label own = mesh_.faceOwner()[facei];
5396 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5398 if (ownLevel < (neiLevel[i]-1))
5402 if (!unrefineCell.test(own))
5408 unrefineCell.unset(own);
5412 else if (neiLevel[i] < (ownLevel-1))
5416 if (unrefineCell.test(own))
5422 unrefineCell.set(own);
5432 Pout<<
"hexRef8::consistentUnrefinement :"
5433 <<
" Changed " << nChanged
5434 <<
" refinement levels due to 2:1 conflicts."
5448 forAll(unrefinePoint, pointi)
5450 if (unrefinePoint.test(pointi))
5452 const labelList& pCells = mesh_.pointCells(pointi);
5456 if (!unrefineCell.test(pCells[j]))
5458 unrefinePoint.unset(pointi);
5470 forAll(unrefinePoint, pointi)
5472 if (unrefinePoint.test(pointi))
5481 forAll(unrefinePoint, pointi)
5483 if (unrefinePoint.test(pointi))
5485 newPointsToUnrefine[nSet++] = pointi;
5489 return newPointsToUnrefine;
5499 if (!history_.active())
5502 <<
"Only call if constructed with history capability"
5508 Pout<<
"hexRef8::setUnrefinement :"
5509 <<
" Checking initial mesh just to make sure" <<
endl;
5513 forAll(cellLevel_, celli)
5515 if (cellLevel_[celli] < 0)
5518 <<
"Illegal cell level " << cellLevel_[celli]
5519 <<
" for cell " << celli
5526 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5529 cellSet cSet(mesh_,
"splitPointCells", splitPointLabels.size());
5531 forAll(splitPointLabels, i)
5533 const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
5535 cSet.insert(pCells);
5539 Pout<<
"hexRef8::setRefinement : Dumping " << pSet.size()
5541 << cSet.size() <<
" cells for unrefinement to" <<
nl
5543 <<
" cellSet " << cSet.objectPath()
5555 forAll(splitPointLabels, i)
5559 splitFaces.insert(
pFaces);
5564 faceRemover_.compatibleRemoves
5572 if (facesToRemove.size() != splitFaces.size())
5576 const_cast<polyMesh&
>(mesh_).setInstance
5578 mesh_.time().timeName()
5581 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5583 faceSet fSet(mesh_,
"splitFaces", splitFaces);
5585 faceSet removeSet(mesh_,
"facesToRemove", facesToRemove);
5590 <<
"Ininitial set of split points to unrefine does not"
5591 <<
" seem to be consistent or not mid points of refined cells"
5592 <<
" splitPoints:" << splitPointLabels.size()
5593 <<
" splitFaces:" << splitFaces.size()
5594 <<
" facesToRemove:" << facesToRemove.size()
5603 forAll(splitPointLabels, i)
5605 label pointi = splitPointLabels[i];
5609 const labelList& pCells = mesh_.pointCells(pointi);
5612 if (pCells.size() != 8)
5615 <<
"splitPoint " << pointi
5616 <<
" should have 8 cells using it. It has " << pCells
5625 label masterCelli =
min(pCells);
5629 label celli = pCells[j];
5631 label region = cellRegion[celli];
5636 <<
"Ininitial set of split points to unrefine does not"
5637 <<
" seem to be consistent or not mid points"
5638 <<
" of refined cells" <<
nl
5639 <<
"cell:" << celli <<
" on splitPoint " << pointi
5640 <<
" has no region to be merged into"
5644 if (masterCelli != cellRegionMaster[region])
5647 <<
"cell:" << celli <<
" on splitPoint:" << pointi
5648 <<
" in region " << region
5649 <<
" has master:" << cellRegionMaster[region]
5650 <<
" which is not the lowest numbered cell"
5651 <<
" among the pointCells:" << pCells
5660 faceRemover_.setRefinement
5671 forAll(splitPointLabels, i)
5673 label pointi = splitPointLabels[i];
5675 const labelList& pCells = mesh_.pointCells(pointi);
5677 label masterCelli =
min(pCells);
5681 cellLevel_[pCells[j]]--;
5684 history_.combineCells(masterCelli, pCells);
5688 setInstance(mesh_.facesInstance());
5698 cellLevel_.write(valid)
5699 && pointLevel_.write(valid)
5700 && level0Edge_.write(valid);
5702 if (history_.active())
5704 writeOk = writeOk && history_.write(valid);
5728 if (
exists(setsDir/
"cellLevel"))
5730 rm(setsDir/
"cellLevel");
5732 if (
exists(setsDir/
"pointLevel"))
5734 rm(setsDir/
"pointLevel");
5736 if (
exists(setsDir/
"level0Edge"))
5738 rm(setsDir/
"level0Edge");