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)];
131 polyTopoChange& meshMod,
144 if ((nei == -1) || (own < nei))
147 newFacei = meshMod.setAction
167 newFacei = meshMod.setAction
171 newFace.reverseFace(),
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);
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)
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],
642 label level = pointLevel_[
f[fp]];
644 if (level < minLevel)
663 label level = pointLevel_[
f[fp]];
665 if (level > maxLevel)
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]];
726 label nexPointi =
f[
f.fcIndex(fp)];
728 str <<
"l " << pointToObjVert[pointi]+1
729 <<
' ' << pointToObjVert[nexPointi]+1 <<
nl;
741 const bool searchForward,
742 const label wantedLevel
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]);
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
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)
1845 if (pointLevel_[pointi] == cellLevel+1)
1847 auto iter = pointFaces.find(pointi);
1852 if (!
pFaces.found(facei))
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()
2044 Foam::hexRef8::hexRef8
2050 const scalar level0Edge
2059 mesh_.facesInstance(),
2072 mesh_.facesInstance(),
2085 mesh_.facesInstance(),
2096 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2103 "refinementHistory",
2104 mesh_.facesInstance(),
2112 faceRemover_(mesh_, GREAT),
2113 savedPointLevel_(0),
2116 if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
2119 <<
"History enabled but number of visible cells in it "
2120 << history_.visibleCells().size()
2121 <<
" is not equal to the number of cells in the mesh "
2127 cellLevel_.size() != mesh_.nCells()
2128 || pointLevel_.size() != mesh_.nPoints()
2132 <<
"Incorrect cellLevel or pointLevel size." <<
endl
2133 <<
"Number of cells in mesh:" << mesh_.nCells()
2134 <<
" does not equal size of cellLevel:" << cellLevel_.size() <<
endl
2135 <<
"Number of points in mesh:" << mesh_.nPoints()
2136 <<
" does not equal size of pointLevel:" << pointLevel_.size()
2141 checkRefinementLevels(-1,
labelList(0));
2154 Foam::hexRef8::hexRef8
2159 const scalar level0Edge
2168 mesh_.facesInstance(),
2181 mesh_.facesInstance(),
2194 mesh_.facesInstance(),
2205 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2212 "refinementHistory",
2213 mesh_.facesInstance(),
2223 faceRemover_(mesh_, GREAT),
2224 savedPointLevel_(0),
2229 cellLevel_.size() != mesh_.nCells()
2230 || pointLevel_.size() != mesh_.nPoints()
2234 <<
"Incorrect cellLevel or pointLevel size." <<
endl
2235 <<
"Number of cells in mesh:" << mesh_.nCells()
2236 <<
" does not equal size of cellLevel:" << cellLevel_.size() <<
endl
2237 <<
"Number of points in mesh:" << mesh_.nPoints()
2238 <<
" does not equal size of pointLevel:" << pointLevel_.size()
2243 checkRefinementLevels(-1,
labelList(0));
2272 label nChanged = faceConsistentRefinement
2283 Pout<<
"hexRef8::consistentRefinement : Changed " << nChanged
2284 <<
" refinement levels due to 2:1 conflicts."
2299 checkWantedRefinementLevels(cellLevel, newCellsToRefine);
2302 return newCellsToRefine;
2314 const label maxFaceDiff,
2317 const label maxPointDiff,
2321 const labelList& faceOwner = mesh_.faceOwner();
2322 const labelList& faceNeighbour = mesh_.faceNeighbour();
2325 if (maxFaceDiff <= 0)
2328 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl
2347 forAll(allCellInfo, celli)
2353 maxFaceDiff*(cellLevel_[celli]+1),
2354 maxFaceDiff*cellLevel_[celli]
2361 label celli = cellsToRefine[i];
2363 allCellInfo[celli].count() = allCellInfo[celli].refinementCount();
2373 int dummyTrackData = 0;
2381 label facei = facesToCheck[i];
2383 if (allFaceInfo[facei].valid(dummyTrackData))
2387 <<
"Argument facesToCheck seems to have duplicate entries!"
2389 <<
"face:" << facei <<
" occurs at positions "
2390 << findIndices(facesToCheck, facei)
2397 if (mesh_.isInternalFace(facei))
2402 const refinementData& neiData = allCellInfo[faceNeighbour[facei]];
2405 label faceRefineCount;
2408 faceCount = neiData.
count() + maxFaceDiff;
2409 faceRefineCount = faceCount + maxFaceDiff;
2413 faceCount = ownData.
count() + maxFaceDiff;
2414 faceRefineCount = faceCount + maxFaceDiff;
2417 seedFaces.append(facei);
2418 seedFacesInfo.append
2426 allFaceInfo[facei] = seedFacesInfo.last();
2430 label faceCount = ownData.
count() + maxFaceDiff;
2431 label faceRefineCount = faceCount + maxFaceDiff;
2433 seedFaces.append(facei);
2434 seedFacesInfo.append
2442 allFaceInfo[facei] = seedFacesInfo.last();
2450 forAll(faceNeighbour, facei)
2453 if (!allFaceInfo[facei].valid(dummyTrackData))
2455 label own = faceOwner[facei];
2456 label nei = faceNeighbour[facei];
2460 if (allCellInfo[own].
count() > allCellInfo[nei].
count())
2462 allFaceInfo[facei].updateFace
2472 seedFacesInfo.append(allFaceInfo[facei]);
2474 else if (allCellInfo[own].
count() < allCellInfo[nei].
count())
2476 allFaceInfo[facei].updateFace
2485 seedFaces.append(facei);
2486 seedFacesInfo.append(allFaceInfo[facei]);
2494 for (
label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
2497 if (!allFaceInfo[facei].valid(dummyTrackData))
2499 label own = faceOwner[facei];
2512 seedFaces.append(facei);
2513 seedFacesInfo.append(faceData);
2531 Pout<<
"hexRef8::consistentSlowRefinement : Seeded "
2532 << seedFaces.size() <<
" faces between cells with different"
2533 <<
" refinement level." <<
endl;
2537 levelCalc.
setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
2539 seedFacesInfo.clear();
2542 levelCalc.
iterate(mesh_.globalData().nTotalFaces()+1);
2550 if (maxPointDiff == -1)
2560 forAll(maxPointCount, pointi)
2562 label& pLevel = maxPointCount[pointi];
2564 const labelList& pCells = mesh_.pointCells(pointi);
2568 pLevel =
max(pLevel, allCellInfo[pCells[i]].
count());
2589 label pointi = pointsToCheck[i];
2594 const labelList& pCells = mesh_.pointCells(pointi);
2598 label celli = pCells[pCelli];
2606 maxPointCount[pointi]
2607 >
cellInfo.count() + maxFaceDiff*maxPointDiff
2615 const cell& cFaces = mesh_.cells()[celli];
2619 label facei = cFaces[cFacei];
2632 if (faceData.
count() > allFaceInfo[facei].count())
2634 changedFacesInfo.insert(facei, faceData);
2641 label nChanged = changedFacesInfo.size();
2651 seedFaces.setCapacity(changedFacesInfo.size());
2652 seedFacesInfo.setCapacity(changedFacesInfo.size());
2656 seedFaces.append(iter.key());
2657 seedFacesInfo.append(iter.val());
2664 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2666 label own = mesh_.faceOwner()[facei];
2669 + (allCellInfo[own].isRefined() ? 1 : 0);
2671 label nei = mesh_.faceNeighbour()[facei];
2674 + (allCellInfo[nei].isRefined() ? 1 : 0);
2676 if (
mag(ownLevel-neiLevel) > 1)
2682 <<
" current level:" << cellLevel_[own]
2683 <<
" current refData:" << allCellInfo[own]
2684 <<
" level after refinement:" << ownLevel
2686 <<
"neighbour cell:" << nei
2687 <<
" current level:" << cellLevel_[nei]
2688 <<
" current refData:" << allCellInfo[nei]
2689 <<
" level after refinement:" << neiLevel
2691 <<
"which does not satisfy 2:1 constraints anymore." <<
nl
2692 <<
"face:" << facei <<
" faceRefData:" << allFaceInfo[facei]
2701 labelList neiLevel(mesh_.nBoundaryFaces());
2702 labelList neiCount(mesh_.nBoundaryFaces());
2703 labelList neiRefCount(mesh_.nBoundaryFaces());
2707 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
2708 neiLevel[i] = cellLevel_[own];
2709 neiCount[i] = allCellInfo[own].count();
2710 neiRefCount[i] = allCellInfo[own].refinementCount();
2721 label facei = i+mesh_.nInternalFaces();
2723 label own = mesh_.faceOwner()[facei];
2726 + (allCellInfo[own].isRefined() ? 1 : 0);
2730 + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
2732 if (
mag(ownLevel - nbrLevel) > 1)
2735 label patchi = mesh_.boundaryMesh().whichPatch(facei);
2738 <<
"Celllevel does not satisfy 2:1 constraint."
2739 <<
" On coupled face "
2741 <<
" refData:" << allFaceInfo[facei]
2742 <<
" on patch " << patchi <<
" "
2743 << mesh_.boundaryMesh()[patchi].name() <<
nl
2744 <<
"owner cell " << own
2745 <<
" current level:" << cellLevel_[own]
2746 <<
" current count:" << allCellInfo[own].count()
2747 <<
" current refCount:"
2748 << allCellInfo[own].refinementCount()
2749 <<
" level after refinement:" << ownLevel
2751 <<
"(coupled) neighbour cell"
2752 <<
" has current level:" << neiLevel[i]
2753 <<
" current count:" << neiCount[i]
2754 <<
" current refCount:" << neiRefCount[i]
2755 <<
" level after refinement:" << nbrLevel
2765 forAll(allCellInfo, celli)
2767 if (allCellInfo[celli].isRefined())
2777 forAll(allCellInfo, celli)
2779 if (allCellInfo[celli].isRefined())
2781 newCellsToRefine[nRefined++] = celli;
2787 Pout<<
"hexRef8::consistentSlowRefinement : From "
2788 << cellsToRefine.size() <<
" to " << newCellsToRefine.size()
2789 <<
" cells to refine." <<
endl;
2792 return newCellsToRefine;
2798 const label maxFaceDiff,
2803 const labelList& faceOwner = mesh_.faceOwner();
2804 const labelList& faceNeighbour = mesh_.faceNeighbour();
2806 if (maxFaceDiff <= 0)
2809 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl
2813 const scalar level0Size = 2*maxFaceDiff*level0EdgeLength();
2828 int dummyTrackData = 0;
2834 label celli = cellsToRefine[i];
2839 mesh_.cellCentres()[celli],
2844 forAll(allCellInfo, celli)
2846 if (!allCellInfo[celli].valid(dummyTrackData))
2851 mesh_.cellCentres()[celli],
2867 label facei = facesToCheck[i];
2869 if (allFaceInfo[facei].valid(dummyTrackData))
2873 <<
"Argument facesToCheck seems to have duplicate entries!"
2875 <<
"face:" << facei <<
" occurs at positions "
2876 << findIndices(facesToCheck, facei)
2880 label own = faceOwner[facei];
2884 allCellInfo[own].valid(dummyTrackData)
2885 ? allCellInfo[own].originLevel()
2889 if (!mesh_.isInternalFace(facei))
2893 const point& fc = mesh_.faceCentres()[facei];
2902 allFaceInfo[facei].updateFace
2914 label nei = faceNeighbour[facei];
2918 allCellInfo[nei].valid(dummyTrackData)
2919 ? allCellInfo[nei].originLevel()
2923 if (ownLevel == neiLevel)
2926 allFaceInfo[facei].updateFace
2935 allFaceInfo[facei].updateFace
2948 allFaceInfo[facei].updateFace
2957 allFaceInfo[facei].updateFace
2969 seedFacesInfo.append(allFaceInfo[facei]);
2976 forAll(faceNeighbour, facei)
2979 if (!allFaceInfo[facei].valid(dummyTrackData))
2981 label own = faceOwner[facei];
2985 allCellInfo[own].valid(dummyTrackData)
2986 ? allCellInfo[own].originLevel()
2990 label nei = faceNeighbour[facei];
2994 allCellInfo[nei].valid(dummyTrackData)
2995 ? allCellInfo[nei].originLevel()
2999 if (ownLevel > neiLevel)
3003 allFaceInfo[facei].updateFace
3012 seedFacesInfo.append(allFaceInfo[facei]);
3014 else if (neiLevel > ownLevel)
3016 seedFaces.append(facei);
3017 allFaceInfo[facei].updateFace
3026 seedFacesInfo.append(allFaceInfo[facei]);
3032 seedFacesInfo.shrink();
3042 mesh_.globalData().nTotalCells()+1,
3083 label celli = cellsToRefine[i];
3085 allCellInfo[celli].originLevel() =
sizeof(
label)*8-2;
3086 allCellInfo[celli].origin() = cc[celli];
3093 forAll(allCellInfo, celli)
3095 label wanted = allCellInfo[celli].wantedLevel(cc[celli]);
3097 if (wanted > cellLevel_[celli]+1)
3102 faceConsistentRefinement(
true, cellLevel_,
refineCell);
3106 label nChanged = faceConsistentRefinement(
true, cellLevel_,
refineCell);
3112 Pout<<
"hexRef8::consistentSlowRefinement2 : Changed " << nChanged
3113 <<
" refinement levels due to 2:1 conflicts."
3128 Pout<<
"hexRef8::consistentSlowRefinement2 : From "
3129 << cellsToRefine.size() <<
" to " << newCellsToRefine.size()
3130 <<
" cells to refine." <<
endl;
3135 cellSet cellsIn(mesh_,
"cellsToRefineIn", cellsToRefine);
3136 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
3137 << cellsIn.
size() <<
" to cellSet "
3142 cellSet cellsOut(mesh_,
"cellsToRefineOut", newCellsToRefine);
3143 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
3144 << cellsOut.
size() <<
" to cellSet "
3153 label nChanged = faceConsistentRefinement(
true, cellLevel_,
refineCell);
3158 mesh_,
"cellsToRefineOut2", newCellsToRefine.size()
3164 cellsOut2.insert(celli);
3167 Pout<<
"hexRef8::consistentSlowRefinement2 : writing "
3168 << cellsOut2.size() <<
" to cellSet "
3169 << cellsOut2.objectPath() <<
endl;
3177 if (
refineCell.test(celli) && !savedRefineCell.test(celli))
3181 <<
"Cell:" << celli <<
" cc:"
3182 << mesh_.cellCentres()[celli]
3183 <<
" was not marked for refinement but does not obey"
3184 <<
" 2:1 constraints."
3191 return newCellsToRefine;
3204 Pout<<
"hexRef8::setRefinement :"
3205 <<
" Checking initial mesh just to make sure" <<
endl;
3214 savedPointLevel_.clear();
3215 savedCellLevel_.clear();
3220 forAll(cellLevel_, celli)
3222 newCellLevel.append(cellLevel_[celli]);
3225 forAll(pointLevel_, pointi)
3227 newPointLevel.append(pointLevel_[pointi]);
3233 Pout<<
"hexRef8::setRefinement :"
3234 <<
" Allocating " << cellLabels.size() <<
" cell midpoints."
3242 labelList cellMidPoint(mesh_.nCells(), -1);
3246 label celli = cellLabels[i];
3248 label anchorPointi = mesh_.faces()[mesh_.cells()[celli][0]][0];
3254 mesh_.cellCentres()[celli],
3261 newPointLevel(cellMidPoint[celli]) = cellLevel_[celli]+1;
3267 cellSet splitCells(mesh_,
"splitCells", cellLabels.size());
3269 forAll(cellMidPoint, celli)
3271 if (cellMidPoint[celli] >= 0)
3273 splitCells.insert(celli);
3277 Pout<<
"hexRef8::setRefinement : Dumping " << splitCells.size()
3278 <<
" cells to split to cellSet " << splitCells.objectPath()
3291 Pout<<
"hexRef8::setRefinement :"
3292 <<
" Allocating edge midpoints."
3301 labelList edgeMidPoint(mesh_.nEdges(), -1);
3304 forAll(cellMidPoint, celli)
3306 if (cellMidPoint[celli] >= 0)
3308 const labelList& cEdges = mesh_.cellEdges(celli);
3312 label edgeI = cEdges[i];
3314 const edge&
e = mesh_.edges()[edgeI];
3318 pointLevel_[
e[0]] <= cellLevel_[celli]
3319 && pointLevel_[
e[1]] <= cellLevel_[celli]
3322 edgeMidPoint[edgeI] = 12345;
3349 forAll(edgeMidPoint, edgeI)
3351 if (edgeMidPoint[edgeI] >= 0)
3354 edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points());
3362 point(-GREAT, -GREAT, -GREAT)
3367 forAll(edgeMidPoint, edgeI)
3369 if (edgeMidPoint[edgeI] >= 0)
3374 const edge&
e = mesh_.edges()[edgeI];
3387 newPointLevel(edgeMidPoint[edgeI]) =
3400 OFstream str(mesh_.time().path()/
"edgeMidPoint.obj");
3402 forAll(edgeMidPoint, edgeI)
3404 if (edgeMidPoint[edgeI] >= 0)
3406 const edge&
e = mesh_.edges()[edgeI];
3412 Pout<<
"hexRef8::setRefinement :"
3413 <<
" Dumping edge centres to split to file " << str.
name() <<
endl;
3423 Pout<<
"hexRef8::setRefinement :"
3424 <<
" Allocating face midpoints."
3430 labelList faceAnchorLevel(mesh_.nFaces());
3432 for (
label facei = 0; facei < mesh_.nFaces(); facei++)
3434 faceAnchorLevel[facei] = faceLevel(facei);
3439 labelList faceMidPoint(mesh_.nFaces(), -1);
3444 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
3446 if (faceAnchorLevel[facei] >= 0)
3448 label own = mesh_.faceOwner()[facei];
3449 label ownLevel = cellLevel_[own];
3450 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3452 label nei = mesh_.faceNeighbour()[facei];
3453 label neiLevel = cellLevel_[nei];
3454 label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
3458 newOwnLevel > faceAnchorLevel[facei]
3459 || newNeiLevel > faceAnchorLevel[facei]
3462 faceMidPoint[facei] = 12345;
3475 labelList newNeiLevel(mesh_.nBoundaryFaces());
3479 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
3480 label ownLevel = cellLevel_[own];
3481 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3483 newNeiLevel[i] = newOwnLevel;
3493 label facei = i+mesh_.nInternalFaces();
3495 if (faceAnchorLevel[facei] >= 0)
3497 label own = mesh_.faceOwner()[facei];
3498 label ownLevel = cellLevel_[own];
3499 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3503 newOwnLevel > faceAnchorLevel[facei]
3504 || newNeiLevel[i] > faceAnchorLevel[facei]
3507 faceMidPoint[facei] = 12345;
3532 mesh_.nBoundaryFaces(),
3533 point(-GREAT, -GREAT, -GREAT)
3538 label facei = i+mesh_.nInternalFaces();
3540 if (faceMidPoint[facei] >= 0)
3542 bFaceMids[i] = mesh_.faceCentres()[facei];
3552 forAll(faceMidPoint, facei)
3554 if (faceMidPoint[facei] >= 0)
3559 const face&
f = mesh_.faces()[facei];
3566 facei < mesh_.nInternalFaces()
3567 ? mesh_.faceCentres()[facei]
3568 : bFaceMids[facei-mesh_.nInternalFaces()]
3578 newPointLevel(faceMidPoint[facei]) = faceAnchorLevel[facei]+1;
3585 faceSet splitFaces(mesh_,
"splitFaces", cellLabels.size());
3587 forAll(faceMidPoint, facei)
3589 if (faceMidPoint[facei] >= 0)
3591 splitFaces.insert(facei);
3595 Pout<<
"hexRef8::setRefinement : Dumping " << splitFaces.size()
3596 <<
" faces to split to faceSet " << splitFaces.objectPath() <<
endl;
3617 Pout<<
"hexRef8::setRefinement :"
3618 <<
" Finding cell anchorPoints (8 per cell)"
3631 forAll(cellMidPoint, celli)
3633 if (cellMidPoint[celli] >= 0)
3635 cellAnchorPoints[celli].setSize(8);
3639 forAll(pointLevel_, pointi)
3641 const labelList& pCells = mesh_.pointCells(pointi);
3645 label celli = pCells[pCelli];
3649 cellMidPoint[celli] >= 0
3650 && pointLevel_[pointi] <= cellLevel_[celli]
3653 if (nAnchorPoints[celli] == 8)
3658 <<
" of level " << cellLevel_[celli]
3659 <<
" uses more than 8 points of equal or"
3660 <<
" lower level" <<
nl
3661 <<
"Points so far:" << cellAnchorPoints[celli]
3665 cellAnchorPoints[celli][nAnchorPoints[celli]++] = pointi;
3670 forAll(cellMidPoint, celli)
3672 if (cellMidPoint[celli] >= 0)
3674 if (nAnchorPoints[celli] != 8)
3678 const labelList& cPoints = mesh_.cellPoints(celli);
3682 <<
" of level " << cellLevel_[celli]
3683 <<
" does not seem to have 8 points of equal or"
3684 <<
" lower level" <<
endl
3685 <<
"cellPoints:" << cPoints <<
endl
3700 Pout<<
"hexRef8::setRefinement :"
3701 <<
" Adding cells (1 per anchorPoint)"
3708 forAll(cellAnchorPoints, celli)
3710 const labelList& cAnchors = cellAnchorPoints[celli];
3712 if (cAnchors.size() == 8)
3714 labelList& cAdded = cellAddedCells[celli];
3721 newCellLevel[celli] = cellLevel_[celli]+1;
3724 for (
label i = 1; i < 8; i++)
3734 mesh_.cellZones().whichZone(celli)
3738 newCellLevel(cAdded[i]) = cellLevel_[celli]+1;
3753 Pout<<
"hexRef8::setRefinement :"
3754 <<
" Marking faces to be handled"
3759 bitSet affectedFace(mesh_.nFaces());
3762 forAll(cellMidPoint, celli)
3764 if (cellMidPoint[celli] >= 0)
3766 const cell& cFaces = mesh_.cells()[celli];
3768 affectedFace.
set(cFaces);
3772 forAll(faceMidPoint, facei)
3774 if (faceMidPoint[facei] >= 0)
3776 affectedFace.set(facei);
3780 forAll(edgeMidPoint, edgeI)
3782 if (edgeMidPoint[edgeI] >= 0)
3784 const labelList& eFaces = mesh_.edgeFaces(edgeI);
3786 affectedFace.
set(eFaces);
3797 Pout<<
"hexRef8::setRefinement : Splitting faces" <<
endl;
3800 forAll(faceMidPoint, facei)
3802 if (faceMidPoint[facei] >= 0 && affectedFace.test(facei))
3808 const face&
f = mesh_.faces()[facei];
3812 bool modifiedFace =
false;
3814 label anchorLevel = faceAnchorLevel[facei];
3822 if (pointLevel_[pointi] <= anchorLevel)
3828 faceVerts.
append(pointi);
3844 faceVerts.
append(faceMidPoint[facei]);
3877 if (mesh_.isInternalFace(facei))
3879 label oldOwn = mesh_.faceOwner()[facei];
3880 label oldNei = mesh_.faceNeighbour()[facei];
3882 checkInternalOrientation
3887 mesh_.cellCentres()[oldOwn],
3888 mesh_.cellCentres()[oldNei],
3894 label oldOwn = mesh_.faceOwner()[facei];
3896 checkBoundaryOrientation
3901 mesh_.cellCentres()[oldOwn],
3902 mesh_.faceCentres()[facei],
3911 modifiedFace =
true;
3913 modFace(meshMod, facei, newFace, own, nei);
3917 addFace(meshMod, facei, newFace, own, nei);
3923 affectedFace.unset(facei);
3933 Pout<<
"hexRef8::setRefinement :"
3934 <<
" Adding edge splits to unsplit faces"
3941 forAll(edgeMidPoint, edgeI)
3943 if (edgeMidPoint[edgeI] >= 0)
3947 const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
3951 label facei = eFaces[i];
3953 if (faceMidPoint[facei] < 0 && affectedFace.test(facei))
3957 const face&
f = mesh_.faces()[facei];
3958 const labelList& fEdges = mesh_.faceEdges
3968 newFaceVerts.append(
f[fp]);
3970 label edgeI = fEdges[fp];
3972 if (edgeMidPoint[edgeI] >= 0)
3974 newFaceVerts.
append(edgeMidPoint[edgeI]);
3983 label anchorFp = findMinLevel(
f);
4000 if (mesh_.isInternalFace(facei))
4002 label oldOwn = mesh_.faceOwner()[facei];
4003 label oldNei = mesh_.faceNeighbour()[facei];
4005 checkInternalOrientation
4010 mesh_.cellCentres()[oldOwn],
4011 mesh_.cellCentres()[oldNei],
4017 label oldOwn = mesh_.faceOwner()[facei];
4019 checkBoundaryOrientation
4024 mesh_.cellCentres()[oldOwn],
4025 mesh_.faceCentres()[facei],
4031 modFace(meshMod, facei, newFace, own, nei);
4034 affectedFace.unset(facei);
4046 Pout<<
"hexRef8::setRefinement :"
4047 <<
" Changing owner/neighbour for otherwise unaffected faces"
4051 forAll(affectedFace, facei)
4053 if (affectedFace.test(facei))
4055 const face&
f = mesh_.faces()[facei];
4059 label anchorFp = findMinLevel(
f);
4073 modFace(meshMod, facei,
f, own, nei);
4076 affectedFace.unset(facei);
4091 Pout<<
"hexRef8::setRefinement :"
4092 <<
" Create new internal faces for split cells"
4096 forAll(cellMidPoint, celli)
4098 if (cellMidPoint[celli] >= 0)
4123 forAll(cellMidPoint, celli)
4125 if (cellMidPoint[celli] >= 0)
4127 minPointi =
min(minPointi, cellMidPoint[celli]);
4128 maxPointi =
max(maxPointi, cellMidPoint[celli]);
4131 forAll(faceMidPoint, facei)
4133 if (faceMidPoint[facei] >= 0)
4135 minPointi =
min(minPointi, faceMidPoint[facei]);
4136 maxPointi =
max(maxPointi, faceMidPoint[facei]);
4139 forAll(edgeMidPoint, edgeI)
4141 if (edgeMidPoint[edgeI] >= 0)
4143 minPointi =
min(minPointi, edgeMidPoint[edgeI]);
4144 maxPointi =
max(maxPointi, edgeMidPoint[edgeI]);
4148 if (minPointi !=
labelMax && minPointi != mesh_.nPoints())
4151 <<
"Added point labels not consecutive to existing mesh points."
4153 <<
"mesh_.nPoints():" << mesh_.nPoints()
4154 <<
" minPointi:" << minPointi
4155 <<
" maxPointi:" << maxPointi
4160 pointLevel_.transfer(newPointLevel);
4161 cellLevel_.transfer(newCellLevel);
4164 setInstance(mesh_.facesInstance());
4171 if (history_.active())
4175 Pout<<
"hexRef8::setRefinement :"
4176 <<
" Updating refinement history to " << cellLevel_.size()
4177 <<
" cells" <<
endl;
4181 history_.resize(cellLevel_.size());
4183 forAll(cellAddedCells, celli)
4185 const labelList& addedCells = cellAddedCells[celli];
4187 if (addedCells.size())
4190 history_.storeSplit(celli, addedCells);
4201 label celli = cellLabels[i];
4203 refinedCells[i].
transfer(cellAddedCells[celli]);
4206 return refinedCells;
4217 savedPointLevel_.clear();
4218 savedPointLevel_.resize(2*pointsToStore.size());
4221 label pointi = pointsToStore[i];
4222 savedPointLevel_.insert(pointi, pointLevel_[pointi]);
4225 savedCellLevel_.
clear();
4226 savedCellLevel_.resize(2*cellsToStore.size());
4229 label celli = cellsToStore[i];
4230 savedCellLevel_.insert(celli, cellLevel_[celli]);
4242 updateMesh(map, dummyMap, dummyMap, dummyMap);
4260 Pout<<
"hexRef8::updateMesh :"
4261 <<
" Updating various lists"
4270 Pout<<
"hexRef8::updateMesh :"
4272 <<
" cellMap:" << map.
cellMap().size()
4273 <<
" nCells:" << mesh_.nCells()
4275 <<
" cellLevel_:" << cellLevel_.size()
4277 <<
" pointMap:" << map.
pointMap().size()
4278 <<
" nPoints:" << mesh_.nPoints()
4280 <<
" pointLevel_:" << pointLevel_.size()
4284 if (reverseCellMap.size() == cellLevel_.size())
4290 reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
4298 forAll(cellMap, newCelli)
4300 label oldCelli = cellMap[newCelli];
4304 newCellLevel[newCelli] = -1;
4308 newCellLevel[newCelli] = cellLevel_[oldCelli];
4318 const label newCelli = iter.key();
4319 const label storedCelli = iter.val();
4321 const auto fnd = savedCellLevel_.cfind(storedCelli);
4326 <<
"Problem : trying to restore old value for new cell "
4327 << newCelli <<
" but cannot find old cell " << storedCelli
4328 <<
" in map of stored values " << savedCellLevel_
4331 cellLevel_[newCelli] = fnd.val();
4352 if (reversePointMap.size() == pointLevel_.size())
4355 reorder(reversePointMap, mesh_.nPoints(), -1, pointLevel_);
4362 labelList newPointLevel(pointMap.size());
4368 if (oldPointi == -1)
4381 newPointLevel[
newPointi] = pointLevel_[oldPointi];
4384 pointLevel_.
transfer(newPointLevel);
4392 const label storedPointi = iter.val();
4394 auto fnd = savedPointLevel_.find(storedPointi);
4399 <<
"Problem : trying to restore old value for new point "
4400 <<
newPointi <<
" but cannot find old point "
4402 <<
" in map of stored values " << savedPointLevel_
4422 if (history_.active())
4424 history_.updateMesh(map);
4428 setInstance(mesh_.facesInstance());
4431 faceRemover_.updateMesh(map);
4434 cellShapesPtr_.clear();
4449 Pout<<
"hexRef8::subset :"
4450 <<
" Updating various lists"
4454 if (history_.active())
4457 <<
"Subsetting will not work in combination with unrefinement."
4459 <<
"Proceed at your own risk." <<
endl;
4467 forAll(cellMap, newCelli)
4469 newCellLevel[newCelli] = cellLevel_[cellMap[newCelli]];
4472 cellLevel_.transfer(newCellLevel);
4474 if (cellLevel_.found(-1))
4478 <<
"cellLevel_ contains illegal value -1 after mapping:"
4486 labelList newPointLevel(pointMap.size());
4493 pointLevel_.transfer(newPointLevel);
4495 if (pointLevel_.found(-1))
4499 <<
"pointLevel_ contains illegal value -1 after mapping:"
4506 if (history_.active())
4508 history_.subset(pointMap,
faceMap, cellMap);
4512 setInstance(mesh_.facesInstance());
4518 cellShapesPtr_.clear();
4527 Pout<<
"hexRef8::distribute :"
4528 <<
" Updating various lists"
4538 if (history_.active())
4540 history_.distribute(map);
4544 faceRemover_.distribute(map);
4547 cellShapesPtr_.clear();
4553 const scalar smallDim = 1
e-6 * mesh_.bounds().mag();
4557 Pout<<
"hexRef8::checkMesh : Using matching tolerance smallDim:"
4558 << smallDim <<
endl;
4571 nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4593 label own = mesh_.faceOwner()[facei];
4594 label bFacei = facei-mesh_.nInternalFaces();
4600 <<
"Faces do not seem to be correct across coupled"
4601 <<
" boundaries" <<
endl
4602 <<
"Coupled face " << facei
4603 <<
" between owner " << own
4604 <<
" on patch " << pp.
name()
4605 <<
" and coupled neighbour " << nei[bFacei]
4606 <<
" has two faces connected to it:"
4625 neiFaceAreas[i] =
mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
4633 label facei = i+mesh_.nInternalFaces();
4635 const scalar magArea =
mag(mesh_.faceAreas()[facei]);
4637 if (
mag(magArea - neiFaceAreas[i]) > smallDim)
4639 const face&
f = mesh_.faces()[facei];
4640 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4642 dumpCell(mesh_.faceOwner()[facei]);
4645 <<
"Faces do not seem to be correct across coupled"
4646 <<
" boundaries" <<
endl
4647 <<
"Coupled face " << facei
4648 <<
" on patch " << patchi
4649 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4651 <<
" has face area:" << magArea
4652 <<
" (coupled) neighbour face area differs:"
4654 <<
" to within tolerance " << smallDim
4664 labelList nVerts(mesh_.nBoundaryFaces());
4668 nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size();
4676 label facei = i+mesh_.nInternalFaces();
4678 const face&
f = mesh_.faces()[facei];
4680 if (
f.size() != nVerts[i])
4682 dumpCell(mesh_.faceOwner()[facei]);
4684 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4687 <<
"Faces do not seem to be correct across coupled"
4688 <<
" boundaries" <<
endl
4689 <<
"Coupled face " << facei
4690 <<
" on patch " << patchi
4691 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4693 <<
" has size:" <<
f.size()
4694 <<
" (coupled) neighbour face has size:"
4706 pointField anchorPoints(mesh_.nBoundaryFaces());
4710 label facei = i+mesh_.nInternalFaces();
4711 const point& fc = mesh_.faceCentres()[facei];
4712 const face&
f = mesh_.faces()[facei];
4713 const vector anchorVec(mesh_.points()[
f[0]] - fc);
4715 anchorPoints[i] = anchorVec;
4725 label facei = i+mesh_.nInternalFaces();
4726 const point& fc = mesh_.faceCentres()[facei];
4727 const face&
f = mesh_.faces()[facei];
4728 const vector anchorVec(mesh_.points()[
f[0]] - fc);
4730 if (
mag(anchorVec - anchorPoints[i]) > smallDim)
4732 dumpCell(mesh_.faceOwner()[facei]);
4734 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4737 <<
"Faces do not seem to be correct across coupled"
4738 <<
" boundaries" <<
endl
4739 <<
"Coupled face " << facei
4740 <<
" on patch " << patchi
4741 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4743 <<
" has anchor vector:" << anchorVec
4744 <<
" (coupled) neighbour face anchor vector differs:"
4746 <<
" to within tolerance " << smallDim
4754 Pout<<
"hexRef8::checkMesh : Returning" <<
endl;
4761 const label maxPointDiff,
4767 Pout<<
"hexRef8::checkRefinementLevels :"
4768 <<
" Checking 2:1 refinement level" <<
endl;
4773 cellLevel_.size() != mesh_.nCells()
4774 || pointLevel_.size() != mesh_.nPoints()
4778 <<
"cellLevel size should be number of cells"
4779 <<
" and pointLevel size should be number of points."<<
nl
4780 <<
"cellLevel:" << cellLevel_.size()
4781 <<
" mesh.nCells():" << mesh_.nCells() <<
nl
4782 <<
"pointLevel:" << pointLevel_.size()
4783 <<
" mesh.nPoints():" << mesh_.nPoints()
4793 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
4795 label own = mesh_.faceOwner()[facei];
4796 label nei = mesh_.faceNeighbour()[facei];
4798 if (
mag(cellLevel_[own] - cellLevel_[nei]) > 1)
4804 <<
"Celllevel does not satisfy 2:1 constraint." <<
nl
4805 <<
"On face " << facei <<
" owner cell " << own
4806 <<
" has refinement " << cellLevel_[own]
4807 <<
" neighbour cell " << nei <<
" has refinement "
4814 labelList neiLevel(mesh_.nBoundaryFaces());
4818 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4820 neiLevel[i] = cellLevel_[own];
4828 label facei = i+mesh_.nInternalFaces();
4830 label own = mesh_.faceOwner()[facei];
4832 if (
mag(cellLevel_[own] - neiLevel[i]) > 1)
4836 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4839 <<
"Celllevel does not satisfy 2:1 constraint."
4840 <<
" On coupled face " << facei
4841 <<
" on patch " << patchi <<
" "
4842 << mesh_.boundaryMesh()[patchi].name()
4843 <<
" owner cell " << own <<
" has refinement "
4845 <<
" (coupled) neighbour cell has refinement "
4868 forAll(syncPointLevel, pointi)
4870 if (pointLevel_[pointi] != syncPointLevel[pointi])
4873 <<
"PointLevel is not consistent across coupled patches."
4875 <<
"point:" << pointi <<
" coord:" << mesh_.points()[pointi]
4876 <<
" has level " << pointLevel_[pointi]
4877 <<
" whereas the coupled point has level "
4878 << syncPointLevel[pointi]
4886 if (maxPointDiff != -1)
4891 forAll(maxPointLevel, pointi)
4893 const labelList& pCells = mesh_.pointCells(pointi);
4895 label& pLevel = maxPointLevel[pointi];
4899 pLevel =
max(pLevel, cellLevel_[pCells[i]]);
4915 label pointi = pointsToCheck[i];
4917 const labelList& pCells = mesh_.pointCells(pointi);
4921 label celli = pCells[i];
4925 mag(cellLevel_[celli]-maxPointLevel[pointi])
4932 <<
"Too big a difference between"
4933 <<
" point-connected cells." <<
nl
4935 <<
" cellLevel:" << cellLevel_[celli]
4936 <<
" uses point:" << pointi
4937 <<
" coord:" << mesh_.points()[pointi]
4938 <<
" which is also used by a cell with level:"
4939 << maxPointLevel[pointi]
5014 if (cellShapesPtr_.empty())
5018 Pout<<
"hexRef8::cellShapes() : calculating splitHex cellShapes."
5019 <<
" cellLevel:" << cellLevel_.size()
5020 <<
" pointLevel:" << pointLevel_.size()
5027 label nSplitHex = 0;
5028 label nUnrecognised = 0;
5030 forAll(cellLevel_, celli)
5032 if (meshShapes[celli].model().index() == 0)
5034 label level = cellLevel_[celli];
5038 bool haveQuads = matchHexShape
5048 cellShapesPtr_()[celli] = degenerateMatcher::match(faces);
5059 Pout<<
"hexRef8::cellShapes() :"
5060 <<
" nCells:" << mesh_.nCells() <<
" of which" <<
nl
5061 <<
" primitive:" << (mesh_.nCells()-nSplitHex-nUnrecognised)
5063 <<
" split-hex:" << nSplitHex <<
nl
5064 <<
" poly :" << nUnrecognised <<
nl
5068 return *cellShapesPtr_;
5076 checkRefinementLevels(-1,
labelList(0));
5081 Pout<<
"hexRef8::getSplitPoints :"
5082 <<
" Calculating unrefineable points" <<
endl;
5086 if (!history_.active())
5089 <<
"Only call if constructed with history capability"
5097 labelList splitMaster(mesh_.nPoints(), -1);
5103 for (
label pointi = 0; pointi < mesh_.nPoints(); pointi++)
5105 const labelList& pCells = mesh_.pointCells(pointi);
5107 if (pCells.size() != 8)
5109 splitMaster[pointi] = -2;
5114 const labelList& visibleCells = history_.visibleCells();
5116 forAll(visibleCells, celli)
5118 const labelList& cPoints = mesh_.cellPoints(celli);
5120 if (visibleCells[celli] != -1 && history_.parentIndex(celli) >= 0)
5122 label parentIndex = history_.parentIndex(celli);
5127 label pointi = cPoints[i];
5129 label masterCelli = splitMaster[pointi];
5131 if (masterCelli == -1)
5138 splitMaster[pointi] = parentIndex;
5139 splitMasterLevel[pointi] = cellLevel_[celli] - 1;
5141 else if (masterCelli == -2)
5147 (masterCelli != parentIndex)
5148 || (splitMasterLevel[pointi] != cellLevel_[celli] - 1)
5153 splitMaster[pointi] = -2;
5162 label pointi = cPoints[i];
5164 splitMaster[pointi] = -2;
5172 label facei = mesh_.nInternalFaces();
5173 facei < mesh_.nFaces();
5177 const face&
f = mesh_.faces()[facei];
5181 splitMaster[
f[fp]] = -2;
5188 label nSplitPoints = 0;
5190 forAll(splitMaster, pointi)
5192 if (splitMaster[pointi] >= 0)
5201 forAll(splitMaster, pointi)
5203 if (splitMaster[pointi] >= 0)
5205 splitPoints[nSplitPoints++] = pointi;
5283 Pout<<
"hexRef8::consistentUnrefinement :"
5284 <<
" Determining 2:1 consistent unrefinement" <<
endl;
5290 <<
"maxSet not implemented yet."
5300 bitSet unrefinePoint(mesh_.nPoints(), pointsToUnrefine);
5307 bitSet unrefineCell(mesh_.nCells());
5309 forAll(unrefinePoint, pointi)
5311 if (unrefinePoint.test(pointi))
5313 const labelList& pCells = mesh_.pointCells(pointi);
5315 unrefineCell.
set(pCells);
5327 for (
label facei = 0; facei < mesh_.nInternalFaces(); facei++)
5329 label own = mesh_.faceOwner()[facei];
5330 label nei = mesh_.faceNeighbour()[facei];
5332 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5333 label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
5335 if (ownLevel < (neiLevel-1))
5342 unrefineCell.set(nei);
5353 if (!unrefineCell.test(own))
5359 unrefineCell.unset(own);
5363 else if (neiLevel < (ownLevel-1))
5367 unrefineCell.set(own);
5371 if (!unrefineCell.test(nei))
5377 unrefineCell.unset(nei);
5385 labelList neiLevel(mesh_.nBoundaryFaces());
5389 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
5391 neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
5399 label facei = i+mesh_.nInternalFaces();
5400 label own = mesh_.faceOwner()[facei];
5401 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5403 if (ownLevel < (neiLevel[i]-1))
5407 if (!unrefineCell.test(own))
5413 unrefineCell.unset(own);
5417 else if (neiLevel[i] < (ownLevel-1))
5421 if (unrefineCell.test(own))
5427 unrefineCell.set(own);
5437 Pout<<
"hexRef8::consistentUnrefinement :"
5438 <<
" Changed " << nChanged
5439 <<
" refinement levels due to 2:1 conflicts."
5453 forAll(unrefinePoint, pointi)
5455 if (unrefinePoint.test(pointi))
5457 const labelList& pCells = mesh_.pointCells(pointi);
5461 if (!unrefineCell.test(pCells[j]))
5463 unrefinePoint.unset(pointi);
5475 forAll(unrefinePoint, pointi)
5477 if (unrefinePoint.test(pointi))
5486 forAll(unrefinePoint, pointi)
5488 if (unrefinePoint.test(pointi))
5490 newPointsToUnrefine[nSet++] = pointi;
5494 return newPointsToUnrefine;
5504 if (!history_.active())
5507 <<
"Only call if constructed with history capability"
5513 Pout<<
"hexRef8::setUnrefinement :"
5514 <<
" Checking initial mesh just to make sure" <<
endl;
5518 forAll(cellLevel_, celli)
5520 if (cellLevel_[celli] < 0)
5523 <<
"Illegal cell level " << cellLevel_[celli]
5524 <<
" for cell " << celli
5531 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5534 cellSet cSet(mesh_,
"splitPointCells", splitPointLabels.size());
5536 forAll(splitPointLabels, i)
5538 const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
5540 cSet.insert(pCells);
5544 Pout<<
"hexRef8::setRefinement : Dumping " << pSet.
size()
5546 << cSet.size() <<
" cells for unrefinement to" <<
nl
5548 <<
" cellSet " << cSet.objectPath()
5560 forAll(splitPointLabels, i)
5564 splitFaces.insert(
pFaces);
5569 faceRemover_.compatibleRemoves
5577 if (facesToRemove.size() != splitFaces.size())
5581 const_cast<polyMesh&
>(mesh_).setInstance
5583 mesh_.time().timeName()
5586 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5588 faceSet fSet(mesh_,
"splitFaces", splitFaces);
5590 faceSet removeSet(mesh_,
"facesToRemove", facesToRemove);
5595 <<
"Ininitial set of split points to unrefine does not"
5596 <<
" seem to be consistent or not mid points of refined cells"
5597 <<
" splitPoints:" << splitPointLabels.size()
5598 <<
" splitFaces:" << splitFaces.size()
5599 <<
" facesToRemove:" << facesToRemove.size()
5608 forAll(splitPointLabels, i)
5610 label pointi = splitPointLabels[i];
5614 const labelList& pCells = mesh_.pointCells(pointi);
5617 if (pCells.size() != 8)
5620 <<
"splitPoint " << pointi
5621 <<
" should have 8 cells using it. It has " << pCells
5634 label celli = pCells[j];
5636 label region = cellRegion[celli];
5641 <<
"Ininitial set of split points to unrefine does not"
5642 <<
" seem to be consistent or not mid points"
5643 <<
" of refined cells" <<
nl
5644 <<
"cell:" << celli <<
" on splitPoint " << pointi
5645 <<
" has no region to be merged into"
5649 if (masterCelli != cellRegionMaster[region])
5652 <<
"cell:" << celli <<
" on splitPoint:" << pointi
5653 <<
" in region " << region
5654 <<
" has master:" << cellRegionMaster[region]
5655 <<
" which is not the lowest numbered cell"
5656 <<
" among the pointCells:" << pCells
5665 faceRemover_.setRefinement
5676 forAll(splitPointLabels, i)
5678 label pointi = splitPointLabels[i];
5680 const labelList& pCells = mesh_.pointCells(pointi);
5686 cellLevel_[pCells[j]]--;
5689 history_.combineCells(masterCelli, pCells);
5693 setInstance(mesh_.facesInstance());
5703 cellLevel_.write(valid)
5704 && pointLevel_.write(valid)
5705 && level0Edge_.write(valid);
5707 if (history_.active())
5709 writeOk = writeOk && history_.write(valid);
5733 if (
exists(setsDir/
"cellLevel"))
5735 rm(setsDir/
"cellLevel");
5737 if (
exists(setsDir/
"pointLevel"))
5739 rm(setsDir/
"pointLevel");
5741 if (
exists(setsDir/
"level0Edge"))
5743 rm(setsDir/
"level0Edge");