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 if (!
pFaces.found(facei))
1881 label facei =
pFaces[pFacei];
1882 const face&
f = mesh_.faces()[facei];
1883 if (mesh_.faceOwner()[facei] == celli)
1885 fourFaces[pFacei] =
f;
1889 fourFaces[pFacei] =
f.reverseFace();
1895 SubList<face>(fourFaces, fourFaces.size()),
1900 if (edgeLoops.size() == 1)
1906 bigFace.meshPoints(),
1907 bigFace.edgeLoops()[0],
1912 if (verts.size() == 4)
1914 quads.append(face(0));
1923 return (quads.size() == 6);
1938 mesh_.facesInstance(),
1951 mesh_.facesInstance(),
1964 mesh_.facesInstance(),
1977 "refinementHistory",
1978 mesh_.facesInstance(),
1985 (readHistory ? mesh_.nCells() : 0)
1987 faceRemover_(mesh_, GREAT),
1988 savedPointLevel_(0),
2003 <<
"History enabled but number of visible cells "
2006 <<
" is not equal to the number of cells in the mesh "
2013 cellLevel_.size() != mesh_.nCells()
2014 || pointLevel_.size() != mesh_.nPoints()
2018 <<
"Restarted from inconsistent cellLevel or pointLevel files."
2022 <<
"Number of cells in mesh:" << mesh_.nCells()
2023 <<
" does not equal size of cellLevel:" << cellLevel_.size() <<
endl
2024 <<
"Number of points in mesh:" << mesh_.nPoints()
2025 <<
" does not equal size of pointLevel:" << pointLevel_.size()
2043 Foam::hexRef8::hexRef8
2049 const scalar level0Edge
2058 mesh_.facesInstance(),
2071 mesh_.facesInstance(),
2084 mesh_.facesInstance(),
2095 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2102 "refinementHistory",
2103 mesh_.facesInstance(),
2111 faceRemover_(mesh_, GREAT),
2112 savedPointLevel_(0),
2115 if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
2118 <<
"History enabled but number of visible cells in it "
2119 << history_.visibleCells().size()
2120 <<
" is not equal to the number of cells in the mesh "
2126 cellLevel_.size() != mesh_.nCells()
2127 || pointLevel_.size() != mesh_.nPoints()
2131 <<
"Incorrect cellLevel or pointLevel size." <<
endl
2132 <<
"Number of cells in mesh:" << mesh_.nCells()
2133 <<
" does not equal size of cellLevel:" << cellLevel_.size() <<
endl
2134 <<
"Number of points in mesh:" << mesh_.nPoints()
2135 <<
" does not equal size of pointLevel:" << pointLevel_.size()
2140 checkRefinementLevels(-1,
labelList(0));
2152 Foam::hexRef8::hexRef8
2157 const scalar level0Edge
2166 mesh_.facesInstance(),
2179 mesh_.facesInstance(),
2192 mesh_.facesInstance(),
2203 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2210 "refinementHistory",
2211 mesh_.facesInstance(),
2221 faceRemover_(mesh_, GREAT),
2222 savedPointLevel_(0),
2227 cellLevel_.size() != mesh_.nCells()
2228 || pointLevel_.size() != mesh_.nPoints()
2232 <<
"Incorrect cellLevel or pointLevel size." <<
endl
2233 <<
"Number of cells in mesh:" << mesh_.nCells()
2234 <<
" does not equal size of cellLevel:" << cellLevel_.size() <<
endl
2235 <<
"Number of points in mesh:" << mesh_.nPoints()
2236 <<
" does not equal size of pointLevel:" << pointLevel_.size()
2241 checkRefinementLevels(-1,
labelList(0));
2270 label nChanged = faceConsistentRefinement
2281 Pout<<
"hexRef8::consistentRefinement : Changed " << nChanged
2282 <<
" refinement levels due to 2:1 conflicts."
2297 checkWantedRefinementLevels(cellLevel, newCellsToRefine);
2300 return newCellsToRefine;
2312 const label maxFaceDiff,
2315 const label maxPointDiff,
2319 const labelList& faceOwner = mesh_.faceOwner();
2320 const labelList& faceNeighbour = mesh_.faceNeighbour();
2323 if (maxFaceDiff <= 0)
2326 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl
2345 forAll(allCellInfo, celli)
2351 maxFaceDiff*(cellLevel_[celli]+1),
2352 maxFaceDiff*cellLevel_[celli]
2359 label celli = cellsToRefine[i];
2361 allCellInfo[celli].count() = allCellInfo[celli].refinementCount();
2371 int dummyTrackData = 0;
2379 label facei = facesToCheck[i];
2381 if (allFaceInfo[facei].valid(dummyTrackData))
2385 <<
"Argument facesToCheck seems to have duplicate entries!"
2387 <<
"face:" << facei <<
" occurs at positions "
2395 if (mesh_.isInternalFace(facei))
2400 const refinementData& neiData = allCellInfo[faceNeighbour[facei]];
2403 label faceRefineCount;
2406 faceCount = neiData.
count() + maxFaceDiff;
2407 faceRefineCount = faceCount + maxFaceDiff;
2411 faceCount = ownData.
count() + maxFaceDiff;
2412 faceRefineCount = faceCount + maxFaceDiff;
2415 seedFaces.append(facei);
2416 seedFacesInfo.append
2424 allFaceInfo[facei] = seedFacesInfo.last();
2428 label faceCount = ownData.
count() + maxFaceDiff;
2429 label faceRefineCount = faceCount + maxFaceDiff;
2431 seedFaces.append(facei);
2432 seedFacesInfo.append
2440 allFaceInfo[facei] = seedFacesInfo.last();
2448 forAll(faceNeighbour, facei)
2451 if (!allFaceInfo[facei].valid(dummyTrackData))
2453 label own = faceOwner[facei];
2454 label nei = faceNeighbour[facei];
2458 if (allCellInfo[own].
count() > allCellInfo[nei].
count())
2460 allFaceInfo[facei].updateFace
2470 seedFacesInfo.append(allFaceInfo[facei]);
2472 else if (allCellInfo[own].
count() < allCellInfo[nei].
count())
2474 allFaceInfo[facei].updateFace
2483 seedFaces.append(facei);
2484 seedFacesInfo.append(allFaceInfo[facei]);
2492 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
2495 if (!allFaceInfo[facei].valid(dummyTrackData))
2497 label own = faceOwner[facei];
2510 seedFaces.append(facei);
2511 seedFacesInfo.append(faceData);
2529 Pout<<
"hexRef8::consistentSlowRefinement : Seeded "
2530 << seedFaces.size() <<
" faces between cells with different"
2531 <<
" refinement level." <<
endl;
2535 levelCalc.
setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
2537 seedFacesInfo.clear();
2540 levelCalc.
iterate(mesh_.globalData().nTotalFaces()+1);
2548 if (maxPointDiff == -1)
2558 forAll(maxPointCount, pointi)
2560 label& pLevel = maxPointCount[pointi];
2562 const labelList& pCells = mesh_.pointCells(pointi);
2566 pLevel =
max(pLevel, allCellInfo[pCells[i]].
count());
2587 label pointi = pointsToCheck[i];
2592 const labelList& pCells = mesh_.pointCells(pointi);
2596 label celli = pCells[pCelli];
2604 maxPointCount[pointi]
2605 >
cellInfo.count() + maxFaceDiff*maxPointDiff
2613 const cell& cFaces = mesh_.cells()[celli];
2617 label facei = cFaces[cFacei];
2630 if (faceData.
count() > allFaceInfo[facei].count())
2632 changedFacesInfo.insert(facei, faceData);
2639 label nChanged = changedFacesInfo.size();
2649 seedFaces.setCapacity(changedFacesInfo.size());
2650 seedFacesInfo.setCapacity(changedFacesInfo.size());
2654 seedFaces.append(iter.key());
2655 seedFacesInfo.append(iter.val());
2662 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2664 label own = mesh_.faceOwner()[facei];
2667 + (allCellInfo[own].isRefined() ? 1 : 0);
2669 label nei = mesh_.faceNeighbour()[facei];
2672 + (allCellInfo[nei].isRefined() ? 1 : 0);
2674 if (
mag(ownLevel-neiLevel) > 1)
2680 <<
" current level:" << cellLevel_[own]
2681 <<
" current refData:" << allCellInfo[own]
2682 <<
" level after refinement:" << ownLevel
2684 <<
"neighbour cell:" << nei
2685 <<
" current level:" << cellLevel_[nei]
2686 <<
" current refData:" << allCellInfo[nei]
2687 <<
" level after refinement:" << neiLevel
2689 <<
"which does not satisfy 2:1 constraints anymore." <<
nl
2690 <<
"face:" << facei <<
" faceRefData:" << allFaceInfo[facei]
2699 labelList neiLevel(mesh_.nBoundaryFaces());
2700 labelList neiCount(mesh_.nBoundaryFaces());
2701 labelList neiRefCount(mesh_.nBoundaryFaces());
2705 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
2706 neiLevel[i] = cellLevel_[own];
2707 neiCount[i] = allCellInfo[own].count();
2708 neiRefCount[i] = allCellInfo[own].refinementCount();
2719 label facei = i+mesh_.nInternalFaces();
2721 label own = mesh_.faceOwner()[facei];
2724 + (allCellInfo[own].isRefined() ? 1 : 0);
2728 + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
2730 if (
mag(ownLevel - nbrLevel) > 1)
2733 label patchi = mesh_.boundaryMesh().whichPatch(facei);
2736 <<
"Celllevel does not satisfy 2:1 constraint."
2737 <<
" On coupled face "
2739 <<
" refData:" << allFaceInfo[facei]
2740 <<
" on patch " << patchi <<
" "
2741 << mesh_.boundaryMesh()[patchi].name() <<
nl
2742 <<
"owner cell " << own
2743 <<
" current level:" << cellLevel_[own]
2744 <<
" current count:" << allCellInfo[own].count()
2745 <<
" current refCount:"
2746 << allCellInfo[own].refinementCount()
2747 <<
" level after refinement:" << ownLevel
2749 <<
"(coupled) neighbour cell"
2750 <<
" has current level:" << neiLevel[i]
2751 <<
" current count:" << neiCount[i]
2752 <<
" current refCount:" << neiRefCount[i]
2753 <<
" level after refinement:" << nbrLevel
2763 forAll(allCellInfo, celli)
2765 if (allCellInfo[celli].isRefined())
2775 forAll(allCellInfo, celli)
2777 if (allCellInfo[celli].isRefined())
2779 newCellsToRefine[nRefined++] = celli;
2785 Pout<<
"hexRef8::consistentSlowRefinement : From "
2786 << cellsToRefine.size() <<
" to " << newCellsToRefine.size()
2787 <<
" cells to refine." <<
endl;
2790 return newCellsToRefine;
2796 const label maxFaceDiff,
2801 const labelList& faceOwner = mesh_.faceOwner();
2802 const labelList& faceNeighbour = mesh_.faceNeighbour();
2804 if (maxFaceDiff <= 0)
2807 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl
2811 const scalar level0Size = 2*maxFaceDiff*level0EdgeLength();
2826 int dummyTrackData = 0;
2832 label celli = cellsToRefine[i];
2837 mesh_.cellCentres()[celli],
2842 forAll(allCellInfo, celli)
2844 if (!allCellInfo[celli].valid(dummyTrackData))
2849 mesh_.cellCentres()[celli],
2865 label facei = facesToCheck[i];
2867 if (allFaceInfo[facei].valid(dummyTrackData))
2871 <<
"Argument facesToCheck seems to have duplicate entries!"
2873 <<
"face:" << facei <<
" occurs at positions "
2878 label own = faceOwner[facei];
2882 allCellInfo[own].valid(dummyTrackData)
2883 ? allCellInfo[own].originLevel()
2887 if (!mesh_.isInternalFace(facei))
2891 const point& fc = mesh_.faceCentres()[facei];
2900 allFaceInfo[facei].updateFace
2912 label nei = faceNeighbour[facei];
2916 allCellInfo[nei].valid(dummyTrackData)
2917 ? allCellInfo[nei].originLevel()
2921 if (ownLevel == neiLevel)
2924 allFaceInfo[facei].updateFace
2933 allFaceInfo[facei].updateFace
2946 allFaceInfo[facei].updateFace
2955 allFaceInfo[facei].updateFace
2967 seedFacesInfo.append(allFaceInfo[facei]);
2974 forAll(faceNeighbour, facei)
2977 if (!allFaceInfo[facei].valid(dummyTrackData))
2979 label own = faceOwner[facei];
2983 allCellInfo[own].valid(dummyTrackData)
2984 ? allCellInfo[own].originLevel()
2988 label nei = faceNeighbour[facei];
2992 allCellInfo[nei].valid(dummyTrackData)
2993 ? allCellInfo[nei].originLevel()
2997 if (ownLevel > neiLevel)
3001 allFaceInfo[facei].updateFace
3010 seedFacesInfo.append(allFaceInfo[facei]);
3012 else if (neiLevel > ownLevel)
3014 seedFaces.append(facei);
3015 allFaceInfo[facei].updateFace
3024 seedFacesInfo.append(allFaceInfo[facei]);
3030 seedFacesInfo.shrink();
3040 mesh_.globalData().nTotalCells()+1,
3081 label celli = cellsToRefine[i];
3083 allCellInfo[celli].originLevel() =
sizeof(label)*8-2;
3084 allCellInfo[celli].origin() = cc[celli];
3091 forAll(allCellInfo, celli)
3093 label wanted = allCellInfo[celli].wantedLevel(cc[celli]);
3095 if (wanted > cellLevel_[celli]+1)
3100 faceConsistentRefinement(
true, cellLevel_,
refineCell);
3104 label nChanged = faceConsistentRefinement(
true, cellLevel_,
refineCell);
3110 Pout<<
"hexRef8::consistentSlowRefinement2 : Changed " << nChanged
3111 <<
" refinement levels due to 2:1 conflicts."
3126 Pout<<
"hexRef8::consistentSlowRefinement2 : From "
3127 << cellsToRefine.size() <<
" to " << newCellsToRefine.size()
3128 <<
" cells to refine." <<
endl;
3133 cellSet cellsIn(mesh_,
"cellsToRefineIn", cellsToRefine);
3134 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
3135 << cellsIn.
size() <<
" to cellSet "
3140 cellSet cellsOut(mesh_,
"cellsToRefineOut", newCellsToRefine);
3141 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
3142 << cellsOut.
size() <<
" to cellSet "
3151 label nChanged = faceConsistentRefinement(
true, cellLevel_,
refineCell);
3156 mesh_,
"cellsToRefineOut2", newCellsToRefine.size()
3162 cellsOut2.insert(celli);
3165 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
3166 << cellsOut2.size() <<
" to cellSet "
3167 << cellsOut2.objectPath() <<
endl;
3175 if (
refineCell.test(celli) && !savedRefineCell.test(celli))
3179 <<
"Cell:" << celli <<
" cc:"
3180 << mesh_.cellCentres()[celli]
3181 <<
" was not marked for refinement but does not obey"
3182 <<
" 2:1 constraints."
3189 return newCellsToRefine;
3202 Pout<<
"hexRef8::setRefinement :"
3203 <<
" Checking initial mesh just to make sure" <<
endl;
3212 savedPointLevel_.clear();
3213 savedCellLevel_.clear();
3218 forAll(cellLevel_, celli)
3220 newCellLevel.append(cellLevel_[celli]);
3223 forAll(pointLevel_, pointi)
3225 newPointLevel.append(pointLevel_[pointi]);
3231 Pout<<
"hexRef8::setRefinement :"
3232 <<
" Allocating " << cellLabels.size() <<
" cell midpoints."
3240 labelList cellMidPoint(mesh_.nCells(), -1);
3244 label celli = cellLabels[i];
3246 label anchorPointi = mesh_.faces()[mesh_.cells()[celli][0]][0];
3252 mesh_.cellCentres()[celli],
3259 newPointLevel(cellMidPoint[celli]) = cellLevel_[celli]+1;
3265 cellSet splitCells(mesh_,
"splitCells", cellLabels.size());
3267 forAll(cellMidPoint, celli)
3269 if (cellMidPoint[celli] >= 0)
3271 splitCells.insert(celli);
3275 Pout<<
"hexRef8::setRefinement : Dumping " << splitCells.size()
3276 <<
" cells to split to cellSet " << splitCells.objectPath()
3289 Pout<<
"hexRef8::setRefinement :"
3290 <<
" Allocating edge midpoints."
3299 labelList edgeMidPoint(mesh_.nEdges(), -1);
3302 forAll(cellMidPoint, celli)
3304 if (cellMidPoint[celli] >= 0)
3306 const labelList& cEdges = mesh_.cellEdges(celli);
3310 label edgeI = cEdges[i];
3312 const edge&
e = mesh_.edges()[edgeI];
3316 pointLevel_[
e[0]] <= cellLevel_[celli]
3317 && pointLevel_[
e[1]] <= cellLevel_[celli]
3320 edgeMidPoint[edgeI] = 12345;
3347 forAll(edgeMidPoint, edgeI)
3349 if (edgeMidPoint[edgeI] >= 0)
3352 edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points());
3360 point(-GREAT, -GREAT, -GREAT)
3365 forAll(edgeMidPoint, edgeI)
3367 if (edgeMidPoint[edgeI] >= 0)
3372 const edge&
e = mesh_.edges()[edgeI];
3385 newPointLevel(edgeMidPoint[edgeI]) =
3398 OFstream str(mesh_.time().path()/
"edgeMidPoint.obj");
3400 forAll(edgeMidPoint, edgeI)
3402 if (edgeMidPoint[edgeI] >= 0)
3404 const edge&
e = mesh_.edges()[edgeI];
3410 Pout<<
"hexRef8::setRefinement :"
3411 <<
" Dumping edge centres to split to file " << str.
name() <<
endl;
3421 Pout<<
"hexRef8::setRefinement :"
3422 <<
" Allocating face midpoints."
3428 labelList faceAnchorLevel(mesh_.nFaces());
3430 for (label facei = 0; facei < mesh_.nFaces(); facei++)
3432 faceAnchorLevel[facei] = faceLevel(facei);
3437 labelList faceMidPoint(mesh_.nFaces(), -1);
3442 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
3444 if (faceAnchorLevel[facei] >= 0)
3446 label own = mesh_.faceOwner()[facei];
3447 label ownLevel = cellLevel_[own];
3448 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3450 label nei = mesh_.faceNeighbour()[facei];
3451 label neiLevel = cellLevel_[nei];
3452 label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
3456 newOwnLevel > faceAnchorLevel[facei]
3457 || newNeiLevel > faceAnchorLevel[facei]
3460 faceMidPoint[facei] = 12345;
3473 labelList newNeiLevel(mesh_.nBoundaryFaces());
3477 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
3478 label ownLevel = cellLevel_[own];
3479 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3481 newNeiLevel[i] = newOwnLevel;
3491 label facei = i+mesh_.nInternalFaces();
3493 if (faceAnchorLevel[facei] >= 0)
3495 label own = mesh_.faceOwner()[facei];
3496 label ownLevel = cellLevel_[own];
3497 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3501 newOwnLevel > faceAnchorLevel[facei]
3502 || newNeiLevel[i] > faceAnchorLevel[facei]
3505 faceMidPoint[facei] = 12345;
3530 mesh_.nBoundaryFaces(),
3531 point(-GREAT, -GREAT, -GREAT)
3536 label facei = i+mesh_.nInternalFaces();
3538 if (faceMidPoint[facei] >= 0)
3540 bFaceMids[i] = mesh_.faceCentres()[facei];
3550 forAll(faceMidPoint, facei)
3552 if (faceMidPoint[facei] >= 0)
3557 const face&
f = mesh_.faces()[facei];
3564 facei < mesh_.nInternalFaces()
3565 ? mesh_.faceCentres()[facei]
3566 : bFaceMids[facei-mesh_.nInternalFaces()]
3576 newPointLevel(faceMidPoint[facei]) = faceAnchorLevel[facei]+1;
3583 faceSet splitFaces(mesh_,
"splitFaces", cellLabels.size());
3585 forAll(faceMidPoint, facei)
3587 if (faceMidPoint[facei] >= 0)
3589 splitFaces.insert(facei);
3593 Pout<<
"hexRef8::setRefinement : Dumping " << splitFaces.size()
3594 <<
" faces to split to faceSet " << splitFaces.objectPath() <<
endl;
3615 Pout<<
"hexRef8::setRefinement :"
3616 <<
" Finding cell anchorPoints (8 per cell)"
3629 forAll(cellMidPoint, celli)
3631 if (cellMidPoint[celli] >= 0)
3633 cellAnchorPoints[celli].setSize(8);
3637 forAll(pointLevel_, pointi)
3639 const labelList& pCells = mesh_.pointCells(pointi);
3643 label celli = pCells[pCelli];
3647 cellMidPoint[celli] >= 0
3648 && pointLevel_[pointi] <= cellLevel_[celli]
3651 if (nAnchorPoints[celli] == 8)
3656 <<
" of level " << cellLevel_[celli]
3657 <<
" uses more than 8 points of equal or"
3658 <<
" lower level" <<
nl
3659 <<
"Points so far:" << cellAnchorPoints[celli]
3663 cellAnchorPoints[celli][nAnchorPoints[celli]++] = pointi;
3668 forAll(cellMidPoint, celli)
3670 if (cellMidPoint[celli] >= 0)
3672 if (nAnchorPoints[celli] != 8)
3676 const labelList& cPoints = mesh_.cellPoints(celli);
3680 <<
" of level " << cellLevel_[celli]
3681 <<
" does not seem to have 8 points of equal or"
3682 <<
" lower level" <<
endl
3683 <<
"cellPoints:" << cPoints <<
endl
3698 Pout<<
"hexRef8::setRefinement :"
3699 <<
" Adding cells (1 per anchorPoint)"
3706 forAll(cellAnchorPoints, celli)
3708 const labelList& cAnchors = cellAnchorPoints[celli];
3710 if (cAnchors.size() == 8)
3712 labelList& cAdded = cellAddedCells[celli];
3719 newCellLevel[celli] = cellLevel_[celli]+1;
3722 for (label i = 1; i < 8; i++)
3732 mesh_.cellZones().whichZone(celli)
3736 newCellLevel(cAdded[i]) = cellLevel_[celli]+1;
3751 Pout<<
"hexRef8::setRefinement :"
3752 <<
" Marking faces to be handled"
3757 bitSet affectedFace(mesh_.nFaces());
3760 forAll(cellMidPoint, celli)
3762 if (cellMidPoint[celli] >= 0)
3764 const cell& cFaces = mesh_.cells()[celli];
3766 affectedFace.
set(cFaces);
3770 forAll(faceMidPoint, facei)
3772 if (faceMidPoint[facei] >= 0)
3774 affectedFace.set(facei);
3778 forAll(edgeMidPoint, edgeI)
3780 if (edgeMidPoint[edgeI] >= 0)
3782 const labelList& eFaces = mesh_.edgeFaces(edgeI);
3784 affectedFace.
set(eFaces);
3795 Pout<<
"hexRef8::setRefinement : Splitting faces" <<
endl;
3798 forAll(faceMidPoint, facei)
3800 if (faceMidPoint[facei] >= 0 && affectedFace.test(facei))
3806 const face&
f = mesh_.faces()[facei];
3810 bool modifiedFace =
false;
3812 label anchorLevel = faceAnchorLevel[facei];
3818 label pointi =
f[fp];
3820 if (pointLevel_[pointi] <= anchorLevel)
3826 faceVerts.
append(pointi);
3842 faceVerts.
append(faceMidPoint[facei]);
3875 if (mesh_.isInternalFace(facei))
3877 label oldOwn = mesh_.faceOwner()[facei];
3878 label oldNei = mesh_.faceNeighbour()[facei];
3880 checkInternalOrientation
3885 mesh_.cellCentres()[oldOwn],
3886 mesh_.cellCentres()[oldNei],
3892 label oldOwn = mesh_.faceOwner()[facei];
3894 checkBoundaryOrientation
3899 mesh_.cellCentres()[oldOwn],
3900 mesh_.faceCentres()[facei],
3909 modifiedFace =
true;
3911 modFace(meshMod, facei, newFace, own, nei);
3915 addFace(meshMod, facei, newFace, own, nei);
3921 affectedFace.unset(facei);
3931 Pout<<
"hexRef8::setRefinement :"
3932 <<
" Adding edge splits to unsplit faces"
3939 forAll(edgeMidPoint, edgeI)
3941 if (edgeMidPoint[edgeI] >= 0)
3945 const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
3949 label facei = eFaces[i];
3951 if (faceMidPoint[facei] < 0 && affectedFace.test(facei))
3955 const face&
f = mesh_.faces()[facei];
3956 const labelList& fEdges = mesh_.faceEdges
3966 newFaceVerts.append(
f[fp]);
3968 label edgeI = fEdges[fp];
3970 if (edgeMidPoint[edgeI] >= 0)
3972 newFaceVerts.
append(edgeMidPoint[edgeI]);
3981 label anchorFp = findMinLevel(
f);
3998 if (mesh_.isInternalFace(facei))
4000 label oldOwn = mesh_.faceOwner()[facei];
4001 label oldNei = mesh_.faceNeighbour()[facei];
4003 checkInternalOrientation
4008 mesh_.cellCentres()[oldOwn],
4009 mesh_.cellCentres()[oldNei],
4015 label oldOwn = mesh_.faceOwner()[facei];
4017 checkBoundaryOrientation
4022 mesh_.cellCentres()[oldOwn],
4023 mesh_.faceCentres()[facei],
4029 modFace(meshMod, facei, newFace, own, nei);
4032 affectedFace.unset(facei);
4044 Pout<<
"hexRef8::setRefinement :"
4045 <<
" Changing owner/neighbour for otherwise unaffected faces"
4049 forAll(affectedFace, facei)
4051 if (affectedFace.test(facei))
4053 const face&
f = mesh_.faces()[facei];
4057 label anchorFp = findMinLevel(
f);
4071 modFace(meshMod, facei,
f, own, nei);
4074 affectedFace.unset(facei);
4089 Pout<<
"hexRef8::setRefinement :"
4090 <<
" Create new internal faces for split cells"
4094 forAll(cellMidPoint, celli)
4096 if (cellMidPoint[celli] >= 0)
4121 forAll(cellMidPoint, celli)
4123 if (cellMidPoint[celli] >= 0)
4125 minPointi =
min(minPointi, cellMidPoint[celli]);
4126 maxPointi =
max(maxPointi, cellMidPoint[celli]);
4129 forAll(faceMidPoint, facei)
4131 if (faceMidPoint[facei] >= 0)
4133 minPointi =
min(minPointi, faceMidPoint[facei]);
4134 maxPointi =
max(maxPointi, faceMidPoint[facei]);
4137 forAll(edgeMidPoint, edgeI)
4139 if (edgeMidPoint[edgeI] >= 0)
4141 minPointi =
min(minPointi, edgeMidPoint[edgeI]);
4142 maxPointi =
max(maxPointi, edgeMidPoint[edgeI]);
4146 if (minPointi !=
labelMax && minPointi != mesh_.nPoints())
4149 <<
"Added point labels not consecutive to existing mesh points."
4151 <<
"mesh_.nPoints():" << mesh_.nPoints()
4152 <<
" minPointi:" << minPointi
4153 <<
" maxPointi:" << maxPointi
4158 pointLevel_.transfer(newPointLevel);
4159 cellLevel_.transfer(newCellLevel);
4162 setInstance(mesh_.facesInstance());
4169 if (history_.active())
4173 Pout<<
"hexRef8::setRefinement :"
4174 <<
" Updating refinement history to " << cellLevel_.size()
4175 <<
" cells" <<
endl;
4179 history_.resize(cellLevel_.size());
4181 forAll(cellAddedCells, celli)
4183 const labelList& addedCells = cellAddedCells[celli];
4185 if (addedCells.size())
4188 history_.storeSplit(celli, addedCells);
4199 label celli = cellLabels[i];
4201 refinedCells[i].
transfer(cellAddedCells[celli]);
4204 return refinedCells;
4215 savedPointLevel_.clear();
4216 savedPointLevel_.resize(2*pointsToStore.size());
4219 label pointi = pointsToStore[i];
4220 savedPointLevel_.insert(pointi, pointLevel_[pointi]);
4223 savedCellLevel_.
clear();
4224 savedCellLevel_.resize(2*cellsToStore.size());
4227 label celli = cellsToStore[i];
4228 savedCellLevel_.insert(celli, cellLevel_[celli]);
4240 updateMesh(map, dummyMap, dummyMap, dummyMap);
4258 Pout<<
"hexRef8::updateMesh :"
4259 <<
" Updating various lists"
4268 Pout<<
"hexRef8::updateMesh :"
4270 <<
" cellMap:" << map.
cellMap().size()
4271 <<
" nCells:" << mesh_.nCells()
4273 <<
" cellLevel_:" << cellLevel_.size()
4275 <<
" pointMap:" << map.
pointMap().size()
4276 <<
" nPoints:" << mesh_.nPoints()
4278 <<
" pointLevel_:" << pointLevel_.size()
4282 if (reverseCellMap.size() == cellLevel_.size())
4288 reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
4296 forAll(cellMap, newCelli)
4298 label oldCelli = cellMap[newCelli];
4302 newCellLevel[newCelli] = -1;
4306 newCellLevel[newCelli] = cellLevel_[oldCelli];
4316 const label newCelli = iter.key();
4317 const label storedCelli = iter.val();
4319 const auto fnd = savedCellLevel_.cfind(storedCelli);
4324 <<
"Problem : trying to restore old value for new cell "
4325 << newCelli <<
" but cannot find old cell " << storedCelli
4326 <<
" in map of stored values " << savedCellLevel_
4329 cellLevel_[newCelli] = fnd.val();
4350 if (reversePointMap.size() == pointLevel_.size())
4353 reorder(reversePointMap, mesh_.nPoints(), -1, pointLevel_);
4360 labelList newPointLevel(pointMap.size());
4364 const label oldPointi = pointMap[
newPointi];
4366 if (oldPointi == -1)
4379 newPointLevel[
newPointi] = pointLevel_[oldPointi];
4382 pointLevel_.
transfer(newPointLevel);
4390 const label storedPointi = iter.val();
4392 auto fnd = savedPointLevel_.find(storedPointi);
4397 <<
"Problem : trying to restore old value for new point "
4398 <<
newPointi <<
" but cannot find old point "
4400 <<
" in map of stored values " << savedPointLevel_
4420 if (history_.active())
4422 history_.updateMesh(map);
4426 setInstance(mesh_.facesInstance());
4429 faceRemover_.updateMesh(map);
4432 cellShapesPtr_.clear();
4447 Pout<<
"hexRef8::subset :"
4448 <<
" Updating various lists"
4452 if (history_.active())
4455 <<
"Subsetting will not work in combination with unrefinement."
4457 <<
"Proceed at your own risk." <<
endl;
4465 forAll(cellMap, newCelli)
4467 newCellLevel[newCelli] = cellLevel_[cellMap[newCelli]];
4470 cellLevel_.transfer(newCellLevel);
4472 if (cellLevel_.found(-1))
4476 <<
"cellLevel_ contains illegal value -1 after mapping:"
4484 labelList newPointLevel(pointMap.size());
4491 pointLevel_.transfer(newPointLevel);
4493 if (pointLevel_.found(-1))
4497 <<
"pointLevel_ contains illegal value -1 after mapping:"
4504 if (history_.active())
4506 history_.subset(pointMap,
faceMap, cellMap);
4510 setInstance(mesh_.facesInstance());
4516 cellShapesPtr_.clear();
4525 Pout<<
"hexRef8::distribute :"
4526 <<
" Updating various lists"
4536 if (history_.active())
4538 history_.distribute(map);
4542 faceRemover_.distribute(map);
4545 cellShapesPtr_.clear();
4551 const scalar smallDim = 1
e-6 * mesh_.bounds().mag();
4555 Pout<<
"hexRef8::checkMesh : Using matching tolerance smallDim:"
4556 << smallDim <<
endl;
4569 nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4587 label facei = pp.
start();
4591 label own = mesh_.faceOwner()[facei];
4592 label bFacei = facei-mesh_.nInternalFaces();
4594 if (!cellToFace.insert(
labelPair(own, nei[bFacei]), facei))
4598 <<
"Faces do not seem to be correct across coupled"
4599 <<
" boundaries" <<
endl
4600 <<
"Coupled face " << facei
4601 <<
" between owner " << own
4602 <<
" on patch " << pp.
name()
4603 <<
" and coupled neighbour " << nei[bFacei]
4604 <<
" has two faces connected to it:"
4606 << cellToFace[
labelPair(own, nei[bFacei])]
4623 neiFaceAreas[i] =
mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
4631 label facei = i+mesh_.nInternalFaces();
4633 const scalar magArea =
mag(mesh_.faceAreas()[facei]);
4635 if (
mag(magArea - neiFaceAreas[i]) > smallDim)
4637 const face&
f = mesh_.faces()[facei];
4638 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4640 dumpCell(mesh_.faceOwner()[facei]);
4643 <<
"Faces do not seem to be correct across coupled"
4644 <<
" boundaries" <<
endl
4645 <<
"Coupled face " << facei
4646 <<
" on patch " << patchi
4647 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4649 <<
" has face area:" << magArea
4650 <<
" (coupled) neighbour face area differs:"
4652 <<
" to within tolerance " << smallDim
4662 labelList nVerts(mesh_.nBoundaryFaces());
4666 nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size();
4674 label facei = i+mesh_.nInternalFaces();
4676 const face&
f = mesh_.faces()[facei];
4678 if (
f.size() != nVerts[i])
4680 dumpCell(mesh_.faceOwner()[facei]);
4682 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4685 <<
"Faces do not seem to be correct across coupled"
4686 <<
" boundaries" <<
endl
4687 <<
"Coupled face " << facei
4688 <<
" on patch " << patchi
4689 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4691 <<
" has size:" <<
f.size()
4692 <<
" (coupled) neighbour face has size:"
4704 pointField anchorPoints(mesh_.nBoundaryFaces());
4708 label facei = i+mesh_.nInternalFaces();
4709 const point& fc = mesh_.faceCentres()[facei];
4710 const face&
f = mesh_.faces()[facei];
4711 const vector anchorVec(mesh_.points()[
f[0]] - fc);
4713 anchorPoints[i] = anchorVec;
4723 label facei = i+mesh_.nInternalFaces();
4724 const point& fc = mesh_.faceCentres()[facei];
4725 const face&
f = mesh_.faces()[facei];
4726 const vector anchorVec(mesh_.points()[
f[0]] - fc);
4728 if (
mag(anchorVec - anchorPoints[i]) > smallDim)
4730 dumpCell(mesh_.faceOwner()[facei]);
4732 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4735 <<
"Faces do not seem to be correct across coupled"
4736 <<
" boundaries" <<
endl
4737 <<
"Coupled face " << facei
4738 <<
" on patch " << patchi
4739 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4741 <<
" has anchor vector:" << anchorVec
4742 <<
" (coupled) neighbour face anchor vector differs:"
4744 <<
" to within tolerance " << smallDim
4752 Pout<<
"hexRef8::checkMesh : Returning" <<
endl;
4759 const label maxPointDiff,
4765 Pout<<
"hexRef8::checkRefinementLevels :"
4766 <<
" Checking 2:1 refinement level" <<
endl;
4771 cellLevel_.size() != mesh_.nCells()
4772 || pointLevel_.size() != mesh_.nPoints()
4776 <<
"cellLevel size should be number of cells"
4777 <<
" and pointLevel size should be number of points."<<
nl
4778 <<
"cellLevel:" << cellLevel_.size()
4779 <<
" mesh.nCells():" << mesh_.nCells() <<
nl
4780 <<
"pointLevel:" << pointLevel_.size()
4781 <<
" mesh.nPoints():" << mesh_.nPoints()
4791 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
4793 label own = mesh_.faceOwner()[facei];
4794 label nei = mesh_.faceNeighbour()[facei];
4796 if (
mag(cellLevel_[own] - cellLevel_[nei]) > 1)
4802 <<
"Celllevel does not satisfy 2:1 constraint." <<
nl
4803 <<
"On face " << facei <<
" owner cell " << own
4804 <<
" has refinement " << cellLevel_[own]
4805 <<
" neighbour cell " << nei <<
" has refinement "
4812 labelList neiLevel(mesh_.nBoundaryFaces());
4816 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4818 neiLevel[i] = cellLevel_[own];
4826 label facei = i+mesh_.nInternalFaces();
4828 label own = mesh_.faceOwner()[facei];
4830 if (
mag(cellLevel_[own] - neiLevel[i]) > 1)
4834 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4837 <<
"Celllevel does not satisfy 2:1 constraint."
4838 <<
" On coupled face " << facei
4839 <<
" on patch " << patchi <<
" "
4840 << mesh_.boundaryMesh()[patchi].name()
4841 <<
" owner cell " << own <<
" has refinement "
4843 <<
" (coupled) neighbour cell has refinement "
4866 forAll(syncPointLevel, pointi)
4868 if (pointLevel_[pointi] != syncPointLevel[pointi])
4871 <<
"PointLevel is not consistent across coupled patches."
4873 <<
"point:" << pointi <<
" coord:" << mesh_.points()[pointi]
4874 <<
" has level " << pointLevel_[pointi]
4875 <<
" whereas the coupled point has level "
4876 << syncPointLevel[pointi]
4884 if (maxPointDiff != -1)
4889 forAll(maxPointLevel, pointi)
4891 const labelList& pCells = mesh_.pointCells(pointi);
4893 label& pLevel = maxPointLevel[pointi];
4897 pLevel =
max(pLevel, cellLevel_[pCells[i]]);
4913 label pointi = pointsToCheck[i];
4915 const labelList& pCells = mesh_.pointCells(pointi);
4919 label celli = pCells[i];
4923 mag(cellLevel_[celli]-maxPointLevel[pointi])
4930 <<
"Too big a difference between"
4931 <<
" point-connected cells." <<
nl
4933 <<
" cellLevel:" << cellLevel_[celli]
4934 <<
" uses point:" << pointi
4935 <<
" coord:" << mesh_.points()[pointi]
4936 <<
" which is also used by a cell with level:"
4937 << maxPointLevel[pointi]
5012 if (!cellShapesPtr_)
5016 Pout<<
"hexRef8::cellShapes() : calculating splitHex cellShapes."
5017 <<
" cellLevel:" << cellLevel_.size()
5018 <<
" pointLevel:" << pointLevel_.size()
5025 label nSplitHex = 0;
5026 label nUnrecognised = 0;
5028 forAll(cellLevel_, celli)
5030 if (meshShapes[celli].model().index() == 0)
5032 label level = cellLevel_[celli];
5036 bool haveQuads = matchHexShape
5046 cellShapesPtr_()[celli] = degenerateMatcher::match(faces);
5057 Pout<<
"hexRef8::cellShapes() :"
5058 <<
" nCells:" << mesh_.nCells() <<
" of which" <<
nl
5059 <<
" primitive:" << (mesh_.nCells()-nSplitHex-nUnrecognised)
5061 <<
" split-hex:" << nSplitHex <<
nl
5062 <<
" poly :" << nUnrecognised <<
nl
5066 return *cellShapesPtr_;
5074 checkRefinementLevels(-1,
labelList(0));
5079 Pout<<
"hexRef8::getSplitPoints :"
5080 <<
" Calculating unrefineable points" <<
endl;
5084 if (!history_.active())
5087 <<
"Only call if constructed with history capability"
5095 labelList splitMaster(mesh_.nPoints(), -1);
5101 for (label pointi = 0; pointi < mesh_.nPoints(); pointi++)
5103 const labelList& pCells = mesh_.pointCells(pointi);
5105 if (pCells.size() != 8)
5107 splitMaster[pointi] = -2;
5112 const labelList& visibleCells = history_.visibleCells();
5114 forAll(visibleCells, celli)
5116 const labelList& cPoints = mesh_.cellPoints(celli);
5118 if (visibleCells[celli] != -1 && history_.parentIndex(celli) >= 0)
5120 label parentIndex = history_.parentIndex(celli);
5125 label pointi = cPoints[i];
5127 label masterCelli = splitMaster[pointi];
5129 if (masterCelli == -1)
5136 splitMaster[pointi] = parentIndex;
5137 splitMasterLevel[pointi] = cellLevel_[celli] - 1;
5139 else if (masterCelli == -2)
5145 (masterCelli != parentIndex)
5146 || (splitMasterLevel[pointi] != cellLevel_[celli] - 1)
5151 splitMaster[pointi] = -2;
5160 label pointi = cPoints[i];
5162 splitMaster[pointi] = -2;
5170 label facei = mesh_.nInternalFaces();
5171 facei < mesh_.nFaces();
5175 const face&
f = mesh_.faces()[facei];
5179 splitMaster[
f[fp]] = -2;
5186 label nSplitPoints = 0;
5188 forAll(splitMaster, pointi)
5190 if (splitMaster[pointi] >= 0)
5199 forAll(splitMaster, pointi)
5201 if (splitMaster[pointi] >= 0)
5203 splitPoints[nSplitPoints++] = pointi;
5281 Pout<<
"hexRef8::consistentUnrefinement :"
5282 <<
" Determining 2:1 consistent unrefinement" <<
endl;
5288 <<
"maxSet not implemented yet."
5298 bitSet unrefinePoint(mesh_.nPoints(), pointsToUnrefine);
5305 bitSet unrefineCell(mesh_.nCells());
5307 forAll(unrefinePoint, pointi)
5309 if (unrefinePoint.test(pointi))
5311 const labelList& pCells = mesh_.pointCells(pointi);
5313 unrefineCell.
set(pCells);
5325 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
5327 label own = mesh_.faceOwner()[facei];
5328 label nei = mesh_.faceNeighbour()[facei];
5330 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5331 label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
5333 if (ownLevel < (neiLevel-1))
5340 unrefineCell.set(nei);
5351 if (!unrefineCell.test(own))
5357 unrefineCell.unset(own);
5361 else if (neiLevel < (ownLevel-1))
5365 unrefineCell.set(own);
5369 if (!unrefineCell.test(nei))
5375 unrefineCell.unset(nei);
5383 labelList neiLevel(mesh_.nBoundaryFaces());
5387 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
5389 neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
5397 label facei = i+mesh_.nInternalFaces();
5398 label own = mesh_.faceOwner()[facei];
5399 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5401 if (ownLevel < (neiLevel[i]-1))
5405 if (!unrefineCell.test(own))
5411 unrefineCell.unset(own);
5415 else if (neiLevel[i] < (ownLevel-1))
5419 if (unrefineCell.test(own))
5425 unrefineCell.set(own);
5435 Pout<<
"hexRef8::consistentUnrefinement :"
5436 <<
" Changed " << nChanged
5437 <<
" refinement levels due to 2:1 conflicts."
5451 forAll(unrefinePoint, pointi)
5453 if (unrefinePoint.test(pointi))
5455 const labelList& pCells = mesh_.pointCells(pointi);
5459 if (!unrefineCell.test(pCells[j]))
5461 unrefinePoint.unset(pointi);
5473 forAll(unrefinePoint, pointi)
5475 if (unrefinePoint.test(pointi))
5484 forAll(unrefinePoint, pointi)
5486 if (unrefinePoint.test(pointi))
5488 newPointsToUnrefine[nSet++] = pointi;
5492 return newPointsToUnrefine;
5502 if (!history_.active())
5505 <<
"Only call if constructed with history capability"
5511 Pout<<
"hexRef8::setUnrefinement :"
5512 <<
" Checking initial mesh just to make sure" <<
endl;
5516 forAll(cellLevel_, celli)
5518 if (cellLevel_[celli] < 0)
5521 <<
"Illegal cell level " << cellLevel_[celli]
5522 <<
" for cell " << celli
5529 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5532 cellSet cSet(mesh_,
"splitPointCells", splitPointLabels.size());
5534 forAll(splitPointLabels, i)
5536 const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
5538 cSet.insert(pCells);
5542 Pout<<
"hexRef8::setRefinement : Dumping " << pSet.
size()
5544 << cSet.size() <<
" cells for unrefinement to" <<
nl
5546 <<
" cellSet " << cSet.objectPath()
5558 forAll(splitPointLabels, i)
5562 splitFaces.insert(
pFaces);
5567 faceRemover_.compatibleRemoves
5575 if (facesToRemove.size() != splitFaces.size())
5579 const_cast<polyMesh&
>(mesh_).setInstance
5581 mesh_.time().timeName()
5584 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5586 faceSet fSet(mesh_,
"splitFaces", splitFaces);
5588 faceSet removeSet(mesh_,
"facesToRemove", facesToRemove);
5593 <<
"Ininitial set of split points to unrefine does not"
5594 <<
" seem to be consistent or not mid points of refined cells"
5595 <<
" splitPoints:" << splitPointLabels.size()
5596 <<
" splitFaces:" << splitFaces.size()
5597 <<
" facesToRemove:" << facesToRemove.size()
5606 forAll(splitPointLabels, i)
5608 label pointi = splitPointLabels[i];
5612 const labelList& pCells = mesh_.pointCells(pointi);
5615 if (pCells.size() != 8)
5618 <<
"splitPoint " << pointi
5619 <<
" should have 8 cells using it. It has " << pCells
5628 label masterCelli =
min(pCells);
5632 label celli = pCells[j];
5634 label region = cellRegion[celli];
5639 <<
"Ininitial set of split points to unrefine does not"
5640 <<
" seem to be consistent or not mid points"
5641 <<
" of refined cells" <<
nl
5642 <<
"cell:" << celli <<
" on splitPoint " << pointi
5643 <<
" has no region to be merged into"
5647 if (masterCelli != cellRegionMaster[region])
5650 <<
"cell:" << celli <<
" on splitPoint:" << pointi
5651 <<
" in region " << region
5652 <<
" has master:" << cellRegionMaster[region]
5653 <<
" which is not the lowest numbered cell"
5654 <<
" among the pointCells:" << pCells
5663 faceRemover_.setRefinement
5674 forAll(splitPointLabels, i)
5676 label pointi = splitPointLabels[i];
5678 const labelList& pCells = mesh_.pointCells(pointi);
5680 label masterCelli =
min(pCells);
5684 cellLevel_[pCells[j]]--;
5687 history_.combineCells(masterCelli, pCells);
5691 setInstance(mesh_.facesInstance());
5701 cellLevel_.write(valid)
5702 && pointLevel_.write(valid)
5703 && level0Edge_.write(valid);
5705 if (history_.active())
5707 writeOk = writeOk && history_.write(valid);
5731 if (
exists(setsDir/
"cellLevel"))
5733 rm(setsDir/
"cellLevel");
5735 if (
exists(setsDir/
"pointLevel"))
5737 rm(setsDir/
"pointLevel");
5739 if (
exists(setsDir/
"level0Edge"))
5741 rm(setsDir/
"level0Edge");