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,
138 label
patchID, zoneID, zoneFlip;
140 getFaceInfo(facei,
patchID, zoneID, zoneFlip);
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,
299 label
patchID, zoneID, zoneFlip;
301 getFaceInfo(facei,
patchID, zoneID, zoneFlip);
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);
443 Pout<<
"hexRef8::getLevel0EdgeLength() :" 444 <<
" After phase1: Edgelengths (squared) per refinementlevel:" 445 << typEdgeLenSqr <<
endl;
459 const label cLevel = cellLevel_[celli];
461 const labelList& cEdges = mesh_.cellEdges(celli);
465 const edge&
e = mesh_.edges()[cEdges[i]];
467 scalar edgeLenSqr =
magSqr(
e.vec(mesh_.points()));
469 maxEdgeLenSqr[cLevel] =
max(maxEdgeLenSqr[cLevel], edgeLenSqr);
477 Pout<<
"hexRef8::getLevel0EdgeLength() :" 478 <<
" Poor Edgelengths (squared) per refinementlevel:" 479 << maxEdgeLenSqr <<
endl;
486 forAll(typEdgeLenSqr, levelI)
488 if (typEdgeLenSqr[levelI] == GREAT2 && maxEdgeLenSqr[levelI] >= 0)
490 typEdgeLenSqr[levelI] = maxEdgeLenSqr[levelI];
496 Pout<<
"hexRef8::getLevel0EdgeLength() :" 497 <<
" Final Edgelengths (squared) per refinementlevel:" 498 << typEdgeLenSqr <<
endl;
502 scalar level0Size = -1;
504 forAll(typEdgeLenSqr, levelI)
506 scalar lenSqr = typEdgeLenSqr[levelI];
514 Pout<<
"hexRef8::getLevel0EdgeLength() :" 515 <<
" For level:" << levelI
517 <<
" with equivalent level0 len:" << level0Size
524 if (level0Size == -1)
536 Foam::label Foam::hexRef8::getAnchorCell
545 if (cellAnchorPoints[celli].size())
547 label index = cellAnchorPoints[celli].find(pointi);
551 return cellAddedCells[celli][index];
558 const face&
f = mesh_.faces()[facei];
562 label index = cellAnchorPoints[celli].find(
f[fp]);
566 return cellAddedCells[celli][index];
572 Perr<<
"cell:" << celli <<
" anchorPoints:" << cellAnchorPoints[celli]
576 <<
"Could not find point " << pointi
577 <<
" in the anchorPoints for cell " << celli <<
endl 578 <<
"Does your original mesh obey the 2:1 constraint and" 579 <<
" did you use consistentRefinement to make your cells to refine" 580 <<
" obey this constraint as well?" 593 void Foam::hexRef8::getFaceNeighbours
609 mesh_.faceOwner()[facei],
614 if (mesh_.isInternalFace(facei))
620 mesh_.faceNeighbour()[facei],
633 Foam::label Foam::hexRef8::findMinLevel(
const labelList&
f)
const 640 label level = pointLevel_[
f[fp]];
642 if (level < minLevel)
654 Foam::label Foam::hexRef8::findMaxLevel(
const labelList&
f)
const 661 label level = pointLevel_[
f[fp]];
663 if (level > maxLevel)
674 Foam::label Foam::hexRef8::countAnchors
677 const label anchorLevel
684 if (pointLevel_[
f[fp]] <= anchorLevel)
693 void Foam::hexRef8::dumpCell(
const label celli)
const 695 OFstream str(mesh_.time().path()/
"cell_" +
Foam::name(celli) +
".obj");
696 Pout<<
"hexRef8 : Dumping cell as obj to " << str.
name() <<
endl;
698 const cell& cFaces = mesh_.cells()[celli];
700 Map<label> pointToObjVert;
705 const face&
f = mesh_.faces()[cFaces[i]];
709 if (pointToObjVert.insert(
f[fp], objVertI))
719 const face&
f = mesh_.faces()[cFaces[i]];
723 label pointi =
f[fp];
726 str <<
"l " << pointToObjVert[pointi]+1
727 <<
' ' << pointToObjVert[nexPointi]+1 <<
nl;
734 Foam::label Foam::hexRef8::findLevel
739 const bool searchForward,
740 const label wantedLevel
747 label pointi =
f[fp];
749 if (pointLevel_[pointi] < wantedLevel)
751 dumpCell(mesh_.faceOwner()[facei]);
752 if (mesh_.isInternalFace(facei))
754 dumpCell(mesh_.faceNeighbour()[facei]);
760 <<
" startFp:" << startFp
761 <<
" wantedLevel:" << wantedLevel
764 else if (pointLevel_[pointi] == wantedLevel)
779 dumpCell(mesh_.faceOwner()[facei]);
780 if (mesh_.isInternalFace(facei))
782 dumpCell(mesh_.faceNeighbour()[facei]);
788 <<
" startFp:" << startFp
789 <<
" wantedLevel:" << wantedLevel
799 const face&
f = mesh_.faces()[facei];
803 return pointLevel_[
f[findMaxLevel(
f)]];
807 label ownLevel = cellLevel_[mesh_.faceOwner()[facei]];
809 if (countAnchors(
f, ownLevel) == 4)
813 else if (countAnchors(
f, ownLevel+1) == 4)
825 void Foam::hexRef8::checkInternalOrientation
827 polyTopoChange& meshMod,
835 face compactFace(
identity(newFace.size()));
836 pointField compactPoints(meshMod.points(), newFace);
838 const vector areaNorm(compactFace.areaNormal(compactPoints));
840 const vector dir(neiPt - ownPt);
842 if ((dir & areaNorm) < 0)
845 <<
"cell:" << celli <<
" old face:" << facei
846 <<
" newFace:" << newFace <<
endl 847 <<
" coords:" << compactPoints
848 <<
" ownPt:" << ownPt
849 <<
" neiPt:" << neiPt
853 const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
855 const scalar
s = (fcToOwn & areaNorm) / (dir & areaNorm);
857 if (s < 0.1 || s > 0.9)
860 <<
"cell:" << celli <<
" old face:" << facei
861 <<
" newFace:" << newFace <<
endl 862 <<
" coords:" << compactPoints
863 <<
" ownPt:" << ownPt
864 <<
" neiPt:" << neiPt
871 void Foam::hexRef8::checkBoundaryOrientation
873 polyTopoChange& meshMod,
877 const point& boundaryPt,
881 face compactFace(
identity(newFace.size()));
882 pointField compactPoints(meshMod.points(), newFace);
884 const vector areaNorm(compactFace.areaNormal(compactPoints));
886 const vector dir(boundaryPt - ownPt);
888 if ((dir & areaNorm) < 0)
891 <<
"cell:" << celli <<
" old face:" << facei
892 <<
" newFace:" << newFace
893 <<
" coords:" << compactPoints
894 <<
" ownPt:" << ownPt
895 <<
" boundaryPt:" << boundaryPt
899 const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
901 const scalar
s = (fcToOwn & dir) /
magSqr(dir);
903 if (s < 0.7 || s > 1.3)
906 <<
"cell:" << celli <<
" old face:" << facei
907 <<
" newFace:" << newFace
908 <<
" coords:" << compactPoints
909 <<
" ownPt:" << ownPt
910 <<
" boundaryPt:" << boundaryPt
920 void Foam::hexRef8::insertEdgeSplit
925 DynamicList<label>& verts
928 if (
p0 < mesh_.nPoints() && p1 < mesh_.nPoints())
932 if (edgeI != -1 && edgeMidPoint[edgeI] != -1)
934 verts.append(edgeMidPoint[edgeI]);
948 Foam::label Foam::hexRef8::storeMidPointInfo
956 const bool faceOrder,
957 const label edgeMidPointi,
958 const label anchorPointi,
959 const label faceMidPointi,
961 Map<edge>& midPointToAnchors,
962 Map<edge>& midPointToFaceMids,
963 polyTopoChange& meshMod
968 bool changed =
false;
969 bool haveTwoAnchors =
false;
971 auto edgeMidFnd = midPointToAnchors.find(edgeMidPointi);
973 if (!edgeMidFnd.good())
975 midPointToAnchors.insert(edgeMidPointi, edge(anchorPointi, -1));
979 edge&
e = edgeMidFnd.val();
981 if (anchorPointi !=
e[0])
990 if (
e[0] != -1 &&
e[1] != -1)
992 haveTwoAnchors =
true;
996 bool haveTwoFaceMids =
false;
998 auto faceMidFnd = midPointToFaceMids.find(edgeMidPointi);
1000 if (!faceMidFnd.good())
1002 midPointToFaceMids.insert(edgeMidPointi, edge(faceMidPointi, -1));
1006 edge&
e = faceMidFnd.val();
1008 if (faceMidPointi !=
e[0])
1012 e[1] = faceMidPointi;
1017 if (
e[0] != -1 &&
e[1] != -1)
1019 haveTwoFaceMids =
true;
1026 if (changed && haveTwoAnchors && haveTwoFaceMids)
1028 const edge& anchors = midPointToAnchors[edgeMidPointi];
1029 const edge& faceMids = midPointToFaceMids[edgeMidPointi];
1031 label otherFaceMidPointi = faceMids.otherVertex(faceMidPointi);
1038 DynamicList<label> newFaceVerts(4);
1039 if (faceOrder == (mesh_.faceOwner()[facei] == celli))
1041 newFaceVerts.append(faceMidPointi);
1052 newFaceVerts.append(edgeMidPointi);
1062 newFaceVerts.append(otherFaceMidPointi);
1063 newFaceVerts.append(cellMidPoint[celli]);
1067 newFaceVerts.append(otherFaceMidPointi);
1077 newFaceVerts.append(edgeMidPointi);
1087 newFaceVerts.append(faceMidPointi);
1088 newFaceVerts.append(cellMidPoint[celli]);
1092 newFace.transfer(newFaceVerts);
1094 label anchorCell0 = getAnchorCell
1102 label anchorCell1 = getAnchorCell
1108 anchors.otherVertex(anchorPointi)
1115 if (anchorCell0 < anchorCell1)
1120 ownPt = mesh_.points()[anchorPointi];
1121 neiPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1130 ownPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1131 neiPt = mesh_.points()[anchorPointi];
1138 if (anchorCell0 < anchorCell1)
1140 ownPt = mesh_.points()[anchorPointi];
1141 neiPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1145 ownPt = mesh_.points()[anchors.otherVertex(anchorPointi)];
1146 neiPt = mesh_.points()[anchorPointi];
1149 checkInternalOrientation
1160 return addInternalFace
1178 void Foam::hexRef8::createInternalFaces
1188 polyTopoChange& meshMod
1194 const cell& cFaces = mesh_.cells()[celli];
1195 const label cLevel = cellLevel_[celli];
1198 Map<edge> midPointToAnchors(24);
1200 Map<edge> midPointToFaceMids(24);
1203 DynamicList<label> storage;
1207 label nFacesAdded = 0;
1211 label facei = cFaces[i];
1213 const face&
f = mesh_.faces()[facei];
1214 const labelList& fEdges = mesh_.faceEdges(facei, storage);
1219 label faceMidPointi = -1;
1221 label nAnchors = countAnchors(
f, cLevel);
1229 label anchorFp = -1;
1233 if (pointLevel_[
f[fp]] <= cLevel)
1241 label edgeMid = findLevel
1249 label faceMid = findLevel
1258 faceMidPointi =
f[faceMid];
1260 else if (nAnchors == 4)
1265 faceMidPointi = faceMidPoint[facei];
1269 dumpCell(mesh_.faceOwner()[facei]);
1270 if (mesh_.isInternalFace(facei))
1272 dumpCell(mesh_.faceNeighbour()[facei]);
1276 <<
"nAnchors:" << nAnchors
1277 <<
" facei:" << facei
1289 label point0 =
f[fp0];
1291 if (pointLevel_[point0] <= cLevel)
1300 label edgeMidPointi = -1;
1304 if (pointLevel_[
f[fp1]] <= cLevel)
1307 label edgeI = fEdges[fp0];
1309 edgeMidPointi = edgeMidPoint[edgeI];
1311 if (edgeMidPointi == -1)
1315 const labelList& cPoints = mesh_.cellPoints(celli);
1318 <<
"cell:" << celli <<
" cLevel:" << cLevel
1319 <<
" cell points:" << cPoints
1322 <<
" face:" << facei
1326 <<
" faceAnchorLevel:" << faceAnchorLevel[facei]
1327 <<
" faceMidPoint:" << faceMidPoint[facei]
1328 <<
" faceMidPointi:" << faceMidPointi
1336 label edgeMid = findLevel(facei,
f, fp1,
true, cLevel+1);
1338 edgeMidPointi =
f[edgeMid];
1341 label newFacei = storeMidPointInfo
1364 if (nFacesAdded == 12)
1377 if (pointLevel_[
f[fpMin1]] <= cLevel)
1380 label edgeI = fEdges[fpMin1];
1382 edgeMidPointi = edgeMidPoint[edgeI];
1384 if (edgeMidPointi == -1)
1388 const labelList& cPoints = mesh_.cellPoints(celli);
1391 <<
"cell:" << celli <<
" cLevel:" << cLevel
1392 <<
" cell points:" << cPoints
1395 <<
" face:" << facei
1399 <<
" faceAnchorLevel:" << faceAnchorLevel[facei]
1400 <<
" faceMidPoint:" << faceMidPoint[facei]
1401 <<
" faceMidPointi:" << faceMidPointi
1409 label edgeMid = findLevel
1418 edgeMidPointi =
f[edgeMid];
1421 newFacei = storeMidPointInfo
1444 if (nFacesAdded == 12)
1452 if (nFacesAdded == 12)
1460 void Foam::hexRef8::walkFaceToMid
1465 const label startFp,
1466 DynamicList<label>& faceVerts
1469 const face&
f = mesh_.faces()[facei];
1470 const labelList& fEdges = mesh_.faceEdges(facei);
1479 if (edgeMidPoint[fEdges[fp]] >= 0)
1481 faceVerts.
append(edgeMidPoint[fEdges[fp]]);
1486 if (pointLevel_[
f[fp]] <= cLevel)
1492 else if (pointLevel_[
f[fp]] == cLevel+1)
1495 faceVerts.append(
f[fp]);
1499 else if (pointLevel_[
f[fp]] == cLevel+2)
1502 faceVerts.append(
f[fp]);
1509 void Foam::hexRef8::walkFaceFromMid
1514 const label startFp,
1515 DynamicList<label>& faceVerts
1518 const face&
f = mesh_.faces()[facei];
1519 const labelList& fEdges = mesh_.faceEdges(facei);
1525 if (pointLevel_[
f[fp]] <= cLevel)
1530 else if (pointLevel_[
f[fp]] == cLevel+1)
1533 faceVerts.append(
f[fp]);
1536 else if (pointLevel_[
f[fp]] == cLevel+2)
1546 if (edgeMidPoint[fEdges[fp]] >= 0)
1548 faceVerts.append(edgeMidPoint[fEdges[fp]]);
1557 faceVerts.append(
f[fp]);
1564 Foam::label Foam::hexRef8::faceConsistentRefinement
1574 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1576 label own = mesh_.faceOwner()[facei];
1577 label ownLevel = cellLevel[own] + refineCell.get(own);
1579 label nei = mesh_.faceNeighbour()[facei];
1580 label neiLevel = cellLevel[nei] + refineCell.get(nei);
1582 if (ownLevel > (neiLevel+1))
1586 refineCell.set(nei);
1590 refineCell.unset(own);
1594 else if (neiLevel > (ownLevel+1))
1598 refineCell.set(own);
1602 refineCell.unset(nei);
1611 labelList neiLevel(mesh_.nBoundaryFaces());
1615 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1617 neiLevel[i] = cellLevel[own] + refineCell.get(own);
1626 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1627 label ownLevel = cellLevel[own] + refineCell.get(own);
1629 if (ownLevel > (neiLevel[i]+1))
1633 refineCell.unset(own);
1637 else if (neiLevel[i] > (ownLevel+1))
1641 refineCell.set(own);
1652 void Foam::hexRef8::checkWantedRefinementLevels
1658 bitSet refineCell(mesh_.nCells(), cellsToRefine);
1660 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1662 label own = mesh_.faceOwner()[facei];
1663 label ownLevel = cellLevel[own] + refineCell.get(own);
1665 label nei = mesh_.faceNeighbour()[facei];
1666 label neiLevel = cellLevel[nei] + refineCell.get(nei);
1668 if (
mag(ownLevel-neiLevel) > 1)
1674 <<
" current level:" << cellLevel[own]
1675 <<
" level after refinement:" << ownLevel
1677 <<
"neighbour cell:" << nei
1678 <<
" current level:" << cellLevel[nei]
1679 <<
" level after refinement:" << neiLevel
1681 <<
"which does not satisfy 2:1 constraints anymore." 1688 labelList neiLevel(mesh_.nBoundaryFaces());
1692 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
1694 neiLevel[i] = cellLevel[own] + refineCell.get(own);
1703 label facei = i + mesh_.nInternalFaces();
1705 label own = mesh_.faceOwner()[facei];
1706 label ownLevel = cellLevel[own] + refineCell.get(own);
1708 if (
mag(ownLevel - neiLevel[i]) > 1)
1710 label patchi = mesh_.boundaryMesh().whichPatch(facei);
1714 <<
"Celllevel does not satisfy 2:1 constraint." 1715 <<
" On coupled face " 1717 <<
" on patch " << patchi <<
" " 1718 << mesh_.boundaryMesh()[patchi].name()
1719 <<
" owner cell " << own
1720 <<
" current level:" << cellLevel[own]
1721 <<
" level after refinement:" << ownLevel
1723 <<
" (coupled) neighbour cell will get refinement " 1736 Pout<<
"hexRef8::setInstance(const fileName& inst) : " 1737 <<
"Resetting file instance to " << inst <<
endl;
1740 cellLevel_.instance() = inst;
1741 pointLevel_.instance() = inst;
1742 level0Edge_.instance() = inst;
1743 history_.instance() = inst;
1747 void Foam::hexRef8::collectLevelPoints
1751 DynamicList<label>&
points 1756 if (pointLevel_[
f[fp]] <= level)
1764 void Foam::hexRef8::collectLevelPoints
1769 DynamicList<label>&
points 1774 label pointi = meshPoints[
f[fp]];
1775 if (pointLevel_[pointi] <= level)
1784 bool Foam::hexRef8::matchHexShape
1787 const label cellLevel,
1788 DynamicList<face>& quads
1791 const cell& cFaces = mesh_.cells()[celli];
1794 DynamicList<label> verts(4);
1802 label facei = cFaces[i];
1803 const face&
f = mesh_.faces()[facei];
1806 collectLevelPoints(
f, cellLevel, verts);
1807 if (verts.size() == 4)
1809 if (mesh_.faceOwner()[facei] != celli)
1813 quads.emplace_back(verts);
1818 if (quads.size() < 6)
1820 Map<labelList> pointFaces(2*cFaces.size());
1824 label facei = cFaces[i];
1825 const face&
f = mesh_.faces()[facei];
1833 collectLevelPoints(
f, cellLevel, verts);
1834 if (verts.size() == 1)
1838 for (
const label pointi :
f)
1840 if (pointLevel_[pointi] == cellLevel+1)
1842 pointFaces(pointi).push_uniq(facei);
1859 label facei =
pFaces[pFacei];
1860 const face&
f = mesh_.faces()[facei];
1861 if (mesh_.faceOwner()[facei] == celli)
1863 fourFaces[pFacei] =
f;
1867 fourFaces[pFacei] =
f.reverseFace();
1873 SubList<face>(fourFaces),
1878 if (edgeLoops.size() == 1)
1884 bigFace.meshPoints(),
1885 bigFace.edgeLoops()[0],
1890 if (verts.size() == 4)
1892 quads.emplace_back(verts);
1899 return (quads.size() == 6);
1906 Foam::hexRef8::hexRef8(
const polyMesh&
mesh,
const bool readHistory)
1914 mesh_.facesInstance(),
1927 mesh_.facesInstance(),
1940 mesh_.facesInstance(),
1953 "refinementHistory",
1954 mesh_.facesInstance(),
1961 (readHistory ? mesh_.nCells() : 0)
1963 faceRemover_(mesh_, GREAT),
1964 savedPointLevel_(0),
1979 <<
"History enabled but number of visible cells " 1982 <<
" is not equal to the number of cells in the mesh " 1989 cellLevel_.
size() != mesh_.nCells()
1990 || pointLevel_.
size() != mesh_.nPoints()
1994 <<
"Restarted from inconsistent cellLevel or pointLevel files." 1998 <<
"Number of cells in mesh:" << mesh_.nCells()
1999 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2000 <<
"Number of points in mesh:" << mesh_.nPoints()
2001 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2019 Foam::hexRef8::hexRef8
2021 const polyMesh&
mesh,
2024 const refinementHistory& history,
2025 const scalar level0Edge
2034 mesh_.facesInstance(),
2035 polyMesh::meshSubDir,
2047 mesh_.facesInstance(),
2048 polyMesh::meshSubDir,
2060 mesh_.facesInstance(),
2061 polyMesh::meshSubDir,
2071 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2078 "refinementHistory",
2079 mesh_.facesInstance(),
2080 polyMesh::meshSubDir,
2087 faceRemover_(mesh_, GREAT),
2088 savedPointLevel_(0),
2094 <<
"History enabled but number of visible cells in it " 2096 <<
" is not equal to the number of cells in the mesh " 2102 cellLevel_.
size() != mesh_.nCells()
2103 || pointLevel_.
size() != mesh_.nPoints()
2107 <<
"Incorrect cellLevel or pointLevel size." <<
endl 2108 <<
"Number of cells in mesh:" << mesh_.nCells()
2109 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2110 <<
"Number of points in mesh:" << mesh_.nPoints()
2111 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2128 Foam::hexRef8::hexRef8
2130 const polyMesh&
mesh,
2133 const scalar level0Edge
2142 mesh_.facesInstance(),
2143 polyMesh::meshSubDir,
2155 mesh_.facesInstance(),
2156 polyMesh::meshSubDir,
2168 mesh_.facesInstance(),
2169 polyMesh::meshSubDir,
2179 (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
2186 "refinementHistory",
2187 mesh_.facesInstance(),
2188 polyMesh::meshSubDir,
2193 List<refinementHistory::splitCell8>(0),
2197 faceRemover_(mesh_, GREAT),
2198 savedPointLevel_(0),
2203 cellLevel_.
size() != mesh_.nCells()
2204 || pointLevel_.
size() != mesh_.nPoints()
2208 <<
"Incorrect cellLevel or pointLevel size." <<
endl 2209 <<
"Number of cells in mesh:" << mesh_.nCells()
2210 <<
" does not equal size of cellLevel:" << cellLevel_.
size() <<
endl 2211 <<
"Number of points in mesh:" << mesh_.nPoints()
2212 <<
" does not equal size of pointLevel:" << pointLevel_.
size()
2246 label nChanged = faceConsistentRefinement
2257 Pout<<
"hexRef8::consistentRefinement : Changed " << nChanged
2258 <<
" refinement levels due to 2:1 conflicts." 2273 checkWantedRefinementLevels(cellLevel, newCellsToRefine);
2276 return newCellsToRefine;
2288 const label maxFaceDiff,
2291 const label maxPointDiff,
2295 const labelList& faceOwner = mesh_.faceOwner();
2296 const labelList& faceNeighbour = mesh_.faceNeighbour();
2299 if (maxFaceDiff <= 0)
2302 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl 2321 forAll(allCellInfo, celli)
2327 maxFaceDiff*(cellLevel_[celli]+1),
2328 maxFaceDiff*cellLevel_[celli]
2335 label celli = cellsToRefine[i];
2337 allCellInfo[celli].count() = allCellInfo[celli].refinementCount();
2347 int dummyTrackData = 0;
2355 label facei = facesToCheck[i];
2357 if (allFaceInfo[facei].valid(dummyTrackData))
2361 <<
"Argument facesToCheck seems to have duplicate entries!" 2363 <<
"face:" << facei <<
" occurs at positions " 2371 label maxDataCount = ownData.
count();
2373 if (mesh_.isInternalFace(facei))
2378 const refinementData& neiData = allCellInfo[faceNeighbour[facei]];
2380 if (maxDataCount < neiData.
count())
2382 maxDataCount = neiData.
count();
2386 label faceCount = maxDataCount + maxFaceDiff;
2387 label faceRefineCount = faceCount + maxFaceDiff;
2389 seedFaces.push_back(facei);
2390 allFaceInfo[facei] =
2391 seedFacesInfo.emplace_back(faceRefineCount, faceCount);
2398 forAll(faceNeighbour, facei)
2401 if (!allFaceInfo[facei].valid(dummyTrackData))
2403 label own = faceOwner[facei];
2404 label nei = faceNeighbour[facei];
2408 if (allCellInfo[own].
count() > allCellInfo[nei].
count())
2410 allFaceInfo[facei].updateFace
2420 seedFacesInfo.append(allFaceInfo[facei]);
2422 else if (allCellInfo[own].
count() < allCellInfo[nei].
count())
2424 allFaceInfo[facei].updateFace
2433 seedFaces.append(facei);
2434 seedFacesInfo.append(allFaceInfo[facei]);
2442 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
2445 if (!allFaceInfo[facei].valid(dummyTrackData))
2447 label own = faceOwner[facei];
2460 seedFaces.append(facei);
2461 seedFacesInfo.append(faceData);
2479 Pout<<
"hexRef8::consistentSlowRefinement : Seeded " 2480 << seedFaces.size() <<
" faces between cells with different" 2481 <<
" refinement level." <<
endl;
2485 levelCalc.
setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink());
2487 seedFacesInfo.clear();
2490 levelCalc.
iterate(mesh_.globalData().nTotalFaces()+1);
2498 if (maxPointDiff == -1)
2508 forAll(maxPointCount, pointi)
2510 label& pLevel = maxPointCount[pointi];
2512 const labelList& pCells = mesh_.pointCells(pointi);
2516 pLevel =
max(pLevel, allCellInfo[pCells[i]].
count());
2537 label pointi = pointsToCheck[i];
2542 const labelList& pCells = mesh_.pointCells(pointi);
2546 label celli = pCells[pCelli];
2554 maxPointCount[pointi]
2555 >
cellInfo.count() + maxFaceDiff*maxPointDiff
2563 const cell& cFaces = mesh_.cells()[celli];
2567 label facei = cFaces[cFacei];
2580 if (faceData.
count() > allFaceInfo[facei].count())
2582 changedFacesInfo.insert(facei, faceData);
2589 label nChanged = changedFacesInfo.size();
2599 seedFaces.setCapacity(changedFacesInfo.size());
2600 seedFacesInfo.setCapacity(changedFacesInfo.size());
2604 seedFaces.append(iter.key());
2605 seedFacesInfo.append(iter.val());
2612 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2614 label own = mesh_.faceOwner()[facei];
2617 + (allCellInfo[own].isRefined() ? 1 : 0);
2619 label nei = mesh_.faceNeighbour()[facei];
2622 + (allCellInfo[nei].isRefined() ? 1 : 0);
2624 if (
mag(ownLevel-neiLevel) > 1)
2630 <<
" current level:" << cellLevel_[own]
2631 <<
" current refData:" << allCellInfo[own]
2632 <<
" level after refinement:" << ownLevel
2634 <<
"neighbour cell:" << nei
2635 <<
" current level:" << cellLevel_[nei]
2636 <<
" current refData:" << allCellInfo[nei]
2637 <<
" level after refinement:" << neiLevel
2639 <<
"which does not satisfy 2:1 constraints anymore." <<
nl 2640 <<
"face:" << facei <<
" faceRefData:" << allFaceInfo[facei]
2649 labelList neiLevel(mesh_.nBoundaryFaces());
2650 labelList neiCount(mesh_.nBoundaryFaces());
2651 labelList neiRefCount(mesh_.nBoundaryFaces());
2655 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
2656 neiLevel[i] = cellLevel_[own];
2657 neiCount[i] = allCellInfo[own].count();
2658 neiRefCount[i] = allCellInfo[own].refinementCount();
2669 label facei = i+mesh_.nInternalFaces();
2671 label own = mesh_.faceOwner()[facei];
2674 + (allCellInfo[own].isRefined() ? 1 : 0);
2678 + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0);
2680 if (
mag(ownLevel - nbrLevel) > 1)
2683 label patchi = mesh_.boundaryMesh().whichPatch(facei);
2686 <<
"Celllevel does not satisfy 2:1 constraint." 2687 <<
" On coupled face " 2689 <<
" refData:" << allFaceInfo[facei]
2690 <<
" on patch " << patchi <<
" " 2691 << mesh_.boundaryMesh()[patchi].name() <<
nl 2692 <<
"owner cell " << own
2693 <<
" current level:" << cellLevel_[own]
2694 <<
" current count:" << allCellInfo[own].count()
2695 <<
" current refCount:" 2696 << allCellInfo[own].refinementCount()
2697 <<
" level after refinement:" << ownLevel
2699 <<
"(coupled) neighbour cell" 2700 <<
" has current level:" << neiLevel[i]
2701 <<
" current count:" << neiCount[i]
2702 <<
" current refCount:" << neiRefCount[i]
2703 <<
" level after refinement:" << nbrLevel
2713 forAll(allCellInfo, celli)
2715 if (allCellInfo[celli].isRefined())
2725 forAll(allCellInfo, celli)
2727 if (allCellInfo[celli].isRefined())
2729 newCellsToRefine[nRefined++] = celli;
2735 Pout<<
"hexRef8::consistentSlowRefinement : From " 2736 << cellsToRefine.
size() <<
" to " << newCellsToRefine.
size()
2737 <<
" cells to refine." <<
endl;
2740 return newCellsToRefine;
2746 const label maxFaceDiff,
2751 const labelList& faceOwner = mesh_.faceOwner();
2752 const labelList& faceNeighbour = mesh_.faceNeighbour();
2754 if (maxFaceDiff <= 0)
2757 <<
"Illegal maxFaceDiff " << maxFaceDiff <<
nl 2761 const scalar level0Size = 2*maxFaceDiff*level0EdgeLength();
2770 List<refinementDistanceData> allCellInfo(mesh_.nCells());
2773 List<refinementDistanceData> allFaceInfo(mesh_.nFaces());
2776 int dummyTrackData = 0;
2782 label celli = cellsToRefine[i];
2784 allCellInfo[celli] = refinementDistanceData
2787 mesh_.cellCentres()[celli],
2792 forAll(allCellInfo, celli)
2794 if (!allCellInfo[celli].valid(dummyTrackData))
2796 allCellInfo[celli] = refinementDistanceData
2799 mesh_.cellCentres()[celli],
2807 DynamicList<label> seedFaces(mesh_.nFaces()/100);
2809 DynamicList<refinementDistanceData> seedFacesInfo(mesh_.nFaces()/100);
2815 label facei = facesToCheck[i];
2817 if (allFaceInfo[facei].valid(dummyTrackData))
2821 <<
"Argument facesToCheck seems to have duplicate entries!" 2823 <<
"face:" << facei <<
" occurs at positions " 2828 label own = faceOwner[facei];
2832 allCellInfo[own].valid(dummyTrackData)
2833 ? allCellInfo[own].originLevel()
2837 if (!mesh_.isInternalFace(facei))
2841 const point& fc = mesh_.faceCentres()[facei];
2843 refinementDistanceData neiData
2850 allFaceInfo[facei].updateFace
2862 label nei = faceNeighbour[facei];
2866 allCellInfo[nei].valid(dummyTrackData)
2867 ? allCellInfo[nei].originLevel()
2871 if (ownLevel == neiLevel)
2874 allFaceInfo[facei].updateFace
2879 refinementDistanceData(level0Size, cc[nei], neiLevel+1),
2883 allFaceInfo[facei].updateFace
2888 refinementDistanceData(level0Size, cc[own], ownLevel+1),
2896 allFaceInfo[facei].updateFace
2901 refinementDistanceData(level0Size, cc[nei], neiLevel),
2905 allFaceInfo[facei].updateFace
2910 refinementDistanceData(level0Size, cc[own], ownLevel),
2917 seedFacesInfo.append(allFaceInfo[facei]);
2924 forAll(faceNeighbour, facei)
2927 if (!allFaceInfo[facei].valid(dummyTrackData))
2929 label own = faceOwner[facei];
2933 allCellInfo[own].valid(dummyTrackData)
2934 ? allCellInfo[own].originLevel()
2938 label nei = faceNeighbour[facei];
2942 allCellInfo[nei].valid(dummyTrackData)
2943 ? allCellInfo[nei].originLevel()
2947 if (ownLevel > neiLevel)
2951 allFaceInfo[facei].updateFace
2956 refinementDistanceData(level0Size, cc[own], ownLevel),
2960 seedFacesInfo.append(allFaceInfo[facei]);
2962 else if (neiLevel > ownLevel)
2964 seedFaces.append(facei);
2965 allFaceInfo[facei].updateFace
2970 refinementDistanceData(level0Size, cc[nei], neiLevel),
2974 seedFacesInfo.append(allFaceInfo[facei]);
2980 seedFacesInfo.shrink();
2983 FaceCellWave<refinementDistanceData, int> levelCalc
2990 mesh_.globalData().nTotalCells()+1,
3031 label celli = cellsToRefine[i];
3033 allCellInfo[celli].originLevel() =
sizeof(label)*8-2;
3034 allCellInfo[celli].origin() = cc[celli];
3040 bitSet refineCell(mesh_.nCells());
3041 forAll(allCellInfo, celli)
3043 label wanted = allCellInfo[celli].wantedLevel(cc[celli]);
3045 if (wanted > cellLevel_[celli]+1)
3047 refineCell.set(celli);
3050 faceConsistentRefinement(
true, cellLevel_, refineCell);
3054 label nChanged = faceConsistentRefinement(
true, cellLevel_, refineCell);
3056 reduce(nChanged, sumOp<label>());
3060 Pout<<
"hexRef8::consistentSlowRefinement2 : Changed " << nChanged
3061 <<
" refinement levels due to 2:1 conflicts." 3072 labelList newCellsToRefine(refineCell.toc());
3076 Pout<<
"hexRef8::consistentSlowRefinement2 : From " 3077 << cellsToRefine.
size() <<
" to " << newCellsToRefine.size()
3078 <<
" cells to refine." <<
endl;
3083 cellSet cellsIn(mesh_,
"cellsToRefineIn", cellsToRefine);
3084 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3085 << cellsIn.size() <<
" to cellSet " 3086 << cellsIn.objectPath() <<
endl;
3090 cellSet cellsOut(mesh_,
"cellsToRefineOut", newCellsToRefine);
3091 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3092 << cellsOut.size() <<
" to cellSet " 3093 << cellsOut.objectPath() <<
endl;
3098 bitSet refineCell(mesh_.nCells(), newCellsToRefine);
3100 const bitSet savedRefineCell(refineCell);
3101 label nChanged = faceConsistentRefinement(
true, cellLevel_, refineCell);
3106 mesh_,
"cellsToRefineOut2", newCellsToRefine.size()
3108 forAll(refineCell, celli)
3110 if (refineCell.test(celli))
3112 cellsOut2.insert(celli);
3115 Pout<<
"hexRef8::consistentSlowRefinement2 : writing " 3116 << cellsOut2.size() <<
" to cellSet " 3117 << cellsOut2.objectPath() <<
endl;
3123 forAll(refineCell, celli)
3125 if (refineCell.test(celli) && !savedRefineCell.test(celli))
3129 <<
"Cell:" << celli <<
" cc:" 3130 << mesh_.cellCentres()[celli]
3131 <<
" was not marked for refinement but does not obey" 3132 <<
" 2:1 constraints." 3139 return newCellsToRefine;
3152 Pout<<
"hexRef8::setRefinement :" 3153 <<
" Checking initial mesh just to make sure" <<
endl;
3162 savedPointLevel_.clear();
3163 savedCellLevel_.clear();
3167 DynamicList<label> newCellLevel(cellLevel_.size());
3168 forAll(cellLevel_, celli)
3170 newCellLevel.append(cellLevel_[celli]);
3172 DynamicList<label> newPointLevel(pointLevel_.size());
3173 forAll(pointLevel_, pointi)
3175 newPointLevel.append(pointLevel_[pointi]);
3181 Pout<<
"hexRef8::setRefinement :" 3182 <<
" Allocating " << cellLabels.
size() <<
" cell midpoints." 3190 labelList cellMidPoint(mesh_.nCells(), -1);
3194 label celli = cellLabels[i];
3196 label anchorPointi = mesh_.faces()[mesh_.cells()[celli][0]][0];
3202 mesh_.cellCentres()[celli],
3209 newPointLevel(cellMidPoint[celli]) = cellLevel_[celli]+1;
3215 cellSet splitCells(mesh_,
"splitCells", cellLabels.
size());
3217 forAll(cellMidPoint, celli)
3219 if (cellMidPoint[celli] >= 0)
3221 splitCells.insert(celli);
3225 Pout<<
"hexRef8::setRefinement : Dumping " << splitCells.size()
3226 <<
" cells to split to cellSet " << splitCells.objectPath()
3239 Pout<<
"hexRef8::setRefinement :" 3240 <<
" Allocating edge midpoints." 3249 labelList edgeMidPoint(mesh_.nEdges(), -1);
3252 forAll(cellMidPoint, celli)
3254 if (cellMidPoint[celli] >= 0)
3256 const labelList& cEdges = mesh_.cellEdges(celli);
3260 label edgeI = cEdges[i];
3262 const edge&
e = mesh_.edges()[edgeI];
3266 pointLevel_[
e[0]] <= cellLevel_[celli]
3267 && pointLevel_[
e[1]] <= cellLevel_[celli]
3270 edgeMidPoint[edgeI] = 12345;
3297 forAll(edgeMidPoint, edgeI)
3299 if (edgeMidPoint[edgeI] >= 0)
3302 edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points());
3310 point(-GREAT, -GREAT, -GREAT)
3315 forAll(edgeMidPoint, edgeI)
3317 if (edgeMidPoint[edgeI] >= 0)
3322 const edge&
e = mesh_.edges()[edgeI];
3335 newPointLevel(edgeMidPoint[edgeI]) =
3348 OFstream str(mesh_.time().path()/
"edgeMidPoint.obj");
3350 forAll(edgeMidPoint, edgeI)
3352 if (edgeMidPoint[edgeI] >= 0)
3354 const edge&
e = mesh_.edges()[edgeI];
3360 Pout<<
"hexRef8::setRefinement :" 3361 <<
" Dumping edge centres to split to file " << str.
name() <<
endl;
3371 Pout<<
"hexRef8::setRefinement :" 3372 <<
" Allocating face midpoints." 3378 labelList faceAnchorLevel(mesh_.nFaces());
3380 for (label facei = 0; facei < mesh_.nFaces(); facei++)
3382 faceAnchorLevel[facei] = faceLevel(facei);
3387 labelList faceMidPoint(mesh_.nFaces(), -1);
3392 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
3394 if (faceAnchorLevel[facei] >= 0)
3396 label own = mesh_.faceOwner()[facei];
3397 label ownLevel = cellLevel_[own];
3398 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3400 label nei = mesh_.faceNeighbour()[facei];
3401 label neiLevel = cellLevel_[nei];
3402 label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0);
3406 newOwnLevel > faceAnchorLevel[facei]
3407 || newNeiLevel > faceAnchorLevel[facei]
3410 faceMidPoint[facei] = 12345;
3423 labelList newNeiLevel(mesh_.nBoundaryFaces());
3427 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
3428 label ownLevel = cellLevel_[own];
3429 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3431 newNeiLevel[i] = newOwnLevel;
3441 label facei = i+mesh_.nInternalFaces();
3443 if (faceAnchorLevel[facei] >= 0)
3445 label own = mesh_.faceOwner()[facei];
3446 label ownLevel = cellLevel_[own];
3447 label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0);
3451 newOwnLevel > faceAnchorLevel[facei]
3452 || newNeiLevel[i] > faceAnchorLevel[facei]
3455 faceMidPoint[facei] = 12345;
3480 mesh_.nBoundaryFaces(),
3481 point(-GREAT, -GREAT, -GREAT)
3486 label facei = i+mesh_.nInternalFaces();
3488 if (faceMidPoint[facei] >= 0)
3490 bFaceMids[i] = mesh_.faceCentres()[facei];
3500 forAll(faceMidPoint, facei)
3502 if (faceMidPoint[facei] >= 0)
3507 const face&
f = mesh_.faces()[facei];
3514 facei < mesh_.nInternalFaces()
3515 ? mesh_.faceCentres()[facei]
3516 : bFaceMids[facei-mesh_.nInternalFaces()]
3526 newPointLevel(faceMidPoint[facei]) = faceAnchorLevel[facei]+1;
3533 faceSet splitFaces(mesh_,
"splitFaces", cellLabels.
size());
3535 forAll(faceMidPoint, facei)
3537 if (faceMidPoint[facei] >= 0)
3539 splitFaces.insert(facei);
3543 Pout<<
"hexRef8::setRefinement : Dumping " << splitFaces.size()
3544 <<
" faces to split to faceSet " << splitFaces.objectPath() <<
endl;
3565 Pout<<
"hexRef8::setRefinement :" 3566 <<
" Finding cell anchorPoints (8 per cell)" 3579 forAll(cellMidPoint, celli)
3581 if (cellMidPoint[celli] >= 0)
3583 cellAnchorPoints[celli].setSize(8);
3587 forAll(pointLevel_, pointi)
3589 const labelList& pCells = mesh_.pointCells(pointi);
3593 label celli = pCells[pCelli];
3597 cellMidPoint[celli] >= 0
3598 && pointLevel_[pointi] <= cellLevel_[celli]
3601 if (nAnchorPoints[celli] == 8)
3606 <<
" of level " << cellLevel_[celli]
3607 <<
" uses more than 8 points of equal or" 3608 <<
" lower level" <<
nl 3609 <<
"Points so far:" << cellAnchorPoints[celli]
3613 cellAnchorPoints[celli][nAnchorPoints[celli]++] = pointi;
3618 forAll(cellMidPoint, celli)
3620 if (cellMidPoint[celli] >= 0)
3622 if (nAnchorPoints[celli] != 8)
3626 const labelList& cPoints = mesh_.cellPoints(celli);
3630 <<
" of level " << cellLevel_[celli]
3631 <<
" does not seem to have 8 points of equal or" 3632 <<
" lower level" <<
endl 3633 <<
"cellPoints:" << cPoints <<
endl 3648 Pout<<
"hexRef8::setRefinement :" 3649 <<
" Adding cells (1 per anchorPoint)" 3656 forAll(cellAnchorPoints, celli)
3658 const labelList& cAnchors = cellAnchorPoints[celli];
3660 if (cAnchors.size() == 8)
3662 labelList& cAdded = cellAddedCells[celli];
3669 newCellLevel[celli] = cellLevel_[celli]+1;
3672 for (label i = 1; i < 8; i++)
3682 mesh_.cellZones().whichZone(celli)
3686 newCellLevel(cAdded[i]) = cellLevel_[celli]+1;
3701 Pout<<
"hexRef8::setRefinement :" 3702 <<
" Marking faces to be handled" 3707 bitSet affectedFace(mesh_.nFaces());
3710 forAll(cellMidPoint, celli)
3712 if (cellMidPoint[celli] >= 0)
3714 const cell& cFaces = mesh_.cells()[celli];
3716 affectedFace.set(cFaces);
3720 forAll(faceMidPoint, facei)
3722 if (faceMidPoint[facei] >= 0)
3724 affectedFace.set(facei);
3728 forAll(edgeMidPoint, edgeI)
3730 if (edgeMidPoint[edgeI] >= 0)
3732 const labelList& eFaces = mesh_.edgeFaces(edgeI);
3734 affectedFace.
set(eFaces);
3745 Pout<<
"hexRef8::setRefinement : Splitting faces" <<
endl;
3748 forAll(faceMidPoint, facei)
3750 if (faceMidPoint[facei] >= 0 && affectedFace.test(facei))
3756 const face&
f = mesh_.faces()[facei];
3760 bool modifiedFace =
false;
3762 label anchorLevel = faceAnchorLevel[facei];
3768 label pointi =
f[fp];
3770 if (pointLevel_[pointi] <= anchorLevel)
3774 DynamicList<label> faceVerts(4);
3776 faceVerts.append(pointi);
3792 faceVerts.append(faceMidPoint[facei]);
3804 newFace.transfer(faceVerts);
3825 if (mesh_.isInternalFace(facei))
3827 label oldOwn = mesh_.faceOwner()[facei];
3828 label oldNei = mesh_.faceNeighbour()[facei];
3830 checkInternalOrientation
3835 mesh_.cellCentres()[oldOwn],
3836 mesh_.cellCentres()[oldNei],
3842 label oldOwn = mesh_.faceOwner()[facei];
3844 checkBoundaryOrientation
3849 mesh_.cellCentres()[oldOwn],
3850 mesh_.faceCentres()[facei],
3859 modifiedFace =
true;
3861 modFace(meshMod, facei, newFace, own, nei);
3865 addFace(meshMod, facei, newFace, own, nei);
3871 affectedFace.unset(facei);
3881 Pout<<
"hexRef8::setRefinement :" 3882 <<
" Adding edge splits to unsplit faces" 3886 DynamicList<label> eFacesStorage;
3887 DynamicList<label> fEdgesStorage;
3889 forAll(edgeMidPoint, edgeI)
3891 if (edgeMidPoint[edgeI] >= 0)
3895 const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
3899 label facei = eFaces[i];
3901 if (faceMidPoint[facei] < 0 && affectedFace.test(facei))
3905 const face&
f = mesh_.faces()[facei];
3906 const labelList& fEdges = mesh_.faceEdges
3912 DynamicList<label> newFaceVerts(
f.
size());
3916 newFaceVerts.append(
f[fp]);
3918 label edgeI = fEdges[fp];
3920 if (edgeMidPoint[edgeI] >= 0)
3922 newFaceVerts.append(edgeMidPoint[edgeI]);
3927 newFace.transfer(newFaceVerts);
3931 label anchorFp = findMinLevel(
f);
3948 if (mesh_.isInternalFace(facei))
3950 label oldOwn = mesh_.faceOwner()[facei];
3951 label oldNei = mesh_.faceNeighbour()[facei];
3953 checkInternalOrientation
3958 mesh_.cellCentres()[oldOwn],
3959 mesh_.cellCentres()[oldNei],
3965 label oldOwn = mesh_.faceOwner()[facei];
3967 checkBoundaryOrientation
3972 mesh_.cellCentres()[oldOwn],
3973 mesh_.faceCentres()[facei],
3979 modFace(meshMod, facei, newFace, own, nei);
3982 affectedFace.unset(facei);
3994 Pout<<
"hexRef8::setRefinement :" 3995 <<
" Changing owner/neighbour for otherwise unaffected faces" 3999 forAll(affectedFace, facei)
4001 if (affectedFace.test(facei))
4003 const face&
f = mesh_.faces()[facei];
4007 label anchorFp = findMinLevel(
f);
4021 modFace(meshMod, facei,
f, own, nei);
4024 affectedFace.unset(facei);
4039 Pout<<
"hexRef8::setRefinement :" 4040 <<
" Create new internal faces for split cells" 4044 forAll(cellMidPoint, celli)
4046 if (cellMidPoint[celli] >= 0)
4071 forAll(cellMidPoint, celli)
4073 if (cellMidPoint[celli] >= 0)
4075 minPointi =
min(minPointi, cellMidPoint[celli]);
4076 maxPointi =
max(maxPointi, cellMidPoint[celli]);
4079 forAll(faceMidPoint, facei)
4081 if (faceMidPoint[facei] >= 0)
4083 minPointi =
min(minPointi, faceMidPoint[facei]);
4084 maxPointi =
max(maxPointi, faceMidPoint[facei]);
4087 forAll(edgeMidPoint, edgeI)
4089 if (edgeMidPoint[edgeI] >= 0)
4091 minPointi =
min(minPointi, edgeMidPoint[edgeI]);
4092 maxPointi =
max(maxPointi, edgeMidPoint[edgeI]);
4096 if (minPointi !=
labelMax && minPointi != mesh_.nPoints())
4099 <<
"Added point labels not consecutive to existing mesh points." 4101 <<
"mesh_.nPoints():" << mesh_.nPoints()
4102 <<
" minPointi:" << minPointi
4103 <<
" maxPointi:" << maxPointi
4108 pointLevel_.transfer(newPointLevel);
4109 cellLevel_.transfer(newCellLevel);
4112 setInstance(mesh_.facesInstance());
4119 if (history_.active())
4123 Pout<<
"hexRef8::setRefinement :" 4124 <<
" Updating refinement history to " << cellLevel_.size()
4125 <<
" cells" <<
endl;
4129 history_.resize(cellLevel_.size());
4131 forAll(cellAddedCells, celli)
4133 const labelList& addedCells = cellAddedCells[celli];
4135 if (addedCells.size())
4138 history_.storeSplit(celli, addedCells);
4149 label celli = cellLabels[i];
4151 refinedCells[i].
transfer(cellAddedCells[celli]);
4154 return refinedCells;
4165 savedPointLevel_.clear();
4166 savedPointLevel_.reserve(pointsToStore.
size());
4169 label pointi = pointsToStore[i];
4170 savedPointLevel_.insert(pointi, pointLevel_[pointi]);
4173 savedCellLevel_.
clear();
4174 savedCellLevel_.reserve(cellsToStore.
size());
4177 label celli = cellsToStore[i];
4178 savedCellLevel_.insert(celli, cellLevel_[celli]);
4190 updateMesh(map, dummyMap, dummyMap, dummyMap);
4208 Pout<<
"hexRef8::updateMesh :" 4209 <<
" Updating various lists" 4218 Pout<<
"hexRef8::updateMesh :" 4221 <<
" nCells:" << mesh_.nCells()
4223 <<
" cellLevel_:" << cellLevel_.size()
4226 <<
" nPoints:" << mesh_.nPoints()
4228 <<
" pointLevel_:" << pointLevel_.size()
4232 if (reverseCellMap.
size() == cellLevel_.size())
4238 reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_);
4246 forAll(cellMap, newCelli)
4248 label oldCelli = cellMap[newCelli];
4252 newCellLevel[newCelli] = -1;
4256 newCellLevel[newCelli] = cellLevel_[oldCelli];
4266 const label newCelli = iter.key();
4267 const label storedCelli = iter.val();
4269 const auto fnd = savedCellLevel_.cfind(storedCelli);
4274 <<
"Problem : trying to restore old value for new cell " 4275 << newCelli <<
" but cannot find old cell " << storedCelli
4276 <<
" in map of stored values " << savedCellLevel_
4279 cellLevel_[newCelli] = fnd.val();
4300 if (reversePointMap.
size() == pointLevel_.size())
4303 reorder(reversePointMap, mesh_.nPoints(), -1, pointLevel_);
4314 const label oldPointi = pointMap[
newPointi];
4316 if (oldPointi == -1)
4329 newPointLevel[
newPointi] = pointLevel_[oldPointi];
4332 pointLevel_.
transfer(newPointLevel);
4340 const label storedPointi = iter.val();
4342 auto fnd = savedPointLevel_.find(storedPointi);
4347 <<
"Problem : trying to restore old value for new point " 4348 <<
newPointi <<
" but cannot find old point " 4350 <<
" in map of stored values " << savedPointLevel_
4370 if (history_.active())
4372 history_.updateMesh(map);
4376 setInstance(mesh_.facesInstance());
4379 faceRemover_.updateMesh(map);
4382 cellShapesPtr_.clear();
4397 Pout<<
"hexRef8::subset :" 4398 <<
" Updating various lists" 4402 if (history_.active())
4405 <<
"Subsetting will not work in combination with unrefinement." 4407 <<
"Proceed at your own risk." <<
endl;
4415 forAll(cellMap, newCelli)
4417 newCellLevel[newCelli] = cellLevel_[cellMap[newCelli]];
4420 cellLevel_.transfer(newCellLevel);
4422 if (cellLevel_.found(-1))
4426 <<
"cellLevel_ contains illegal value -1 after mapping:" 4441 pointLevel_.transfer(newPointLevel);
4443 if (pointLevel_.found(-1))
4447 <<
"pointLevel_ contains illegal value -1 after mapping:" 4454 if (history_.active())
4456 history_.subset(pointMap,
faceMap, cellMap);
4460 setInstance(mesh_.facesInstance());
4466 cellShapesPtr_.clear();
4475 Pout<<
"hexRef8::distribute :" 4476 <<
" Updating various lists" 4486 if (history_.active())
4488 history_.distribute(map);
4492 faceRemover_.distribute(map);
4495 cellShapesPtr_.clear();
4501 const scalar smallDim = 1
e-6 * mesh_.bounds().mag();
4505 Pout<<
"hexRef8::checkMesh : Using matching tolerance smallDim:" 4506 << smallDim <<
endl;
4519 nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4525 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
4537 label facei =
pp.start();
4541 label own = mesh_.faceOwner()[facei];
4542 label bFacei = facei-mesh_.nInternalFaces();
4544 if (!cellToFace.insert(
labelPair(own, nei[bFacei]), facei))
4548 <<
"Faces do not seem to be correct across coupled" 4549 <<
" boundaries" <<
endl 4550 <<
"Coupled face " << facei
4551 <<
" between owner " << own
4552 <<
" on patch " <<
pp.name()
4553 <<
" and coupled neighbour " << nei[bFacei]
4554 <<
" has two faces connected to it:" 4556 << cellToFace[
labelPair(own, nei[bFacei])]
4573 neiFaceAreas[i] =
mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
4581 label facei = i+mesh_.nInternalFaces();
4583 const scalar magArea =
mag(mesh_.faceAreas()[facei]);
4585 if (
mag(magArea - neiFaceAreas[i]) > smallDim)
4587 const face&
f = mesh_.faces()[facei];
4588 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4590 dumpCell(mesh_.faceOwner()[facei]);
4593 <<
"Faces do not seem to be correct across coupled" 4594 <<
" boundaries" <<
endl 4595 <<
"Coupled face " << facei
4596 <<
" on patch " << patchi
4597 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4598 <<
" coords:" << UIndirectList<point>(mesh_.points(),
f)
4599 <<
" has face area:" << magArea
4600 <<
" (coupled) neighbour face area differs:" 4602 <<
" to within tolerance " << smallDim
4612 labelList nVerts(mesh_.nBoundaryFaces());
4616 nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size();
4624 label facei = i+mesh_.nInternalFaces();
4626 const face&
f = mesh_.faces()[facei];
4628 if (
f.
size() != nVerts[i])
4630 dumpCell(mesh_.faceOwner()[facei]);
4632 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4635 <<
"Faces do not seem to be correct across coupled" 4636 <<
" boundaries" <<
endl 4637 <<
"Coupled face " << facei
4638 <<
" on patch " << patchi
4639 <<
" " << mesh_.boundaryMesh()[patchi].
name()
4640 <<
" coords:" << UIndirectList<point>(mesh_.points(),
f)
4641 <<
" has size:" <<
f.
size()
4642 <<
" (coupled) neighbour face has size:" 4654 pointField anchorPoints(mesh_.nBoundaryFaces());
4658 label facei = i+mesh_.nInternalFaces();
4659 const point& fc = mesh_.faceCentres()[facei];
4660 const face&
f = mesh_.faces()[facei];
4661 const vector anchorVec(mesh_.points()[
f[0]] - fc);
4663 anchorPoints[i] = anchorVec;
4673 label facei = i+mesh_.nInternalFaces();
4674 const point& fc = mesh_.faceCentres()[facei];
4675 const face&
f = mesh_.faces()[facei];
4676 const vector anchorVec(mesh_.points()[
f[0]] - fc);
4678 if (
mag(anchorVec - anchorPoints[i]) > smallDim)
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()
4690 <<
" coords:" << UIndirectList<point>(mesh_.points(),
f)
4691 <<
" has anchor vector:" << anchorVec
4692 <<
" (coupled) neighbour face anchor vector differs:" 4694 <<
" to within tolerance " << smallDim
4702 Pout<<
"hexRef8::checkMesh : Returning" <<
endl;
4709 const label maxPointDiff,
4715 Pout<<
"hexRef8::checkRefinementLevels :" 4716 <<
" Checking 2:1 refinement level" <<
endl;
4721 cellLevel_.size() != mesh_.nCells()
4722 || pointLevel_.size() != mesh_.nPoints()
4726 <<
"cellLevel size should be number of cells" 4727 <<
" and pointLevel size should be number of points."<<
nl 4728 <<
"cellLevel:" << cellLevel_.size()
4729 <<
" mesh.nCells():" << mesh_.nCells() <<
nl 4730 <<
"pointLevel:" << pointLevel_.size()
4731 <<
" mesh.nPoints():" << mesh_.nPoints()
4741 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
4743 label own = mesh_.faceOwner()[facei];
4744 label nei = mesh_.faceNeighbour()[facei];
4746 if (
mag(cellLevel_[own] - cellLevel_[nei]) > 1)
4752 <<
"Celllevel does not satisfy 2:1 constraint." <<
nl 4753 <<
"On face " << facei <<
" owner cell " << own
4754 <<
" has refinement " << cellLevel_[own]
4755 <<
" neighbour cell " << nei <<
" has refinement " 4762 labelList neiLevel(mesh_.nBoundaryFaces());
4766 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
4768 neiLevel[i] = cellLevel_[own];
4776 label facei = i+mesh_.nInternalFaces();
4778 label own = mesh_.faceOwner()[facei];
4780 if (
mag(cellLevel_[own] - neiLevel[i]) > 1)
4784 label patchi = mesh_.boundaryMesh().whichPatch(facei);
4787 <<
"Celllevel does not satisfy 2:1 constraint." 4788 <<
" On coupled face " << facei
4789 <<
" on patch " << patchi <<
" " 4790 << mesh_.boundaryMesh()[patchi].name()
4791 <<
" owner cell " << own <<
" has refinement " 4793 <<
" (coupled) neighbour cell has refinement " 4816 forAll(syncPointLevel, pointi)
4818 if (pointLevel_[pointi] != syncPointLevel[pointi])
4821 <<
"PointLevel is not consistent across coupled patches." 4823 <<
"point:" << pointi <<
" coord:" << mesh_.points()[pointi]
4824 <<
" has level " << pointLevel_[pointi]
4825 <<
" whereas the coupled point has level " 4826 << syncPointLevel[pointi]
4834 if (maxPointDiff != -1)
4839 forAll(maxPointLevel, pointi)
4841 const labelList& pCells = mesh_.pointCells(pointi);
4843 label& pLevel = maxPointLevel[pointi];
4847 pLevel =
max(pLevel, cellLevel_[pCells[i]]);
4863 label pointi = pointsToCheck[i];
4865 const labelList& pCells = mesh_.pointCells(pointi);
4869 label celli = pCells[i];
4873 mag(cellLevel_[celli]-maxPointLevel[pointi])
4880 <<
"Too big a difference between" 4881 <<
" point-connected cells." <<
nl 4883 <<
" cellLevel:" << cellLevel_[celli]
4884 <<
" uses point:" << pointi
4885 <<
" coord:" << mesh_.points()[pointi]
4886 <<
" which is also used by a cell with level:" 4887 << maxPointLevel[pointi]
4962 if (!cellShapesPtr_)
4966 Pout<<
"hexRef8::cellShapes() : calculating splitHex cellShapes." 4967 <<
" cellLevel:" << cellLevel_.size()
4968 <<
" pointLevel:" << pointLevel_.size()
4975 label nSplitHex = 0;
4976 label nUnrecognised = 0;
4978 forAll(cellLevel_, celli)
4980 if (meshShapes[celli].model().index() == 0)
4982 label level = cellLevel_[celli];
4986 bool haveQuads = matchHexShape
4996 cellShapesPtr_()[celli] = degenerateMatcher::match(faces);
5007 Pout<<
"hexRef8::cellShapes() :" 5008 <<
" nCells:" << mesh_.nCells() <<
" of which" <<
nl 5009 <<
" primitive:" << (mesh_.nCells()-nSplitHex-nUnrecognised)
5011 <<
" split-hex:" << nSplitHex <<
nl 5012 <<
" poly :" << nUnrecognised <<
nl 5016 return *cellShapesPtr_;
5024 checkRefinementLevels(-1,
labelList(0));
5029 Pout<<
"hexRef8::getSplitPoints :" 5030 <<
" Calculating unrefineable points" <<
endl;
5034 if (!history_.active())
5037 <<
"Only call if constructed with history capability" 5045 labelList splitMaster(mesh_.nPoints(), -1);
5051 for (label pointi = 0; pointi < mesh_.nPoints(); pointi++)
5053 const labelList& pCells = mesh_.pointCells(pointi);
5055 if (pCells.
size() != 8)
5057 splitMaster[pointi] = -2;
5062 const labelList& visibleCells = history_.visibleCells();
5064 forAll(visibleCells, celli)
5066 const labelList& cPoints = mesh_.cellPoints(celli);
5068 if (visibleCells[celli] != -1 && history_.parentIndex(celli) >= 0)
5070 label parentIndex = history_.parentIndex(celli);
5075 label pointi = cPoints[i];
5077 label masterCelli = splitMaster[pointi];
5079 if (masterCelli == -1)
5086 splitMaster[pointi] = parentIndex;
5087 splitMasterLevel[pointi] = cellLevel_[celli] - 1;
5089 else if (masterCelli == -2)
5095 (masterCelli != parentIndex)
5096 || (splitMasterLevel[pointi] != cellLevel_[celli] - 1)
5101 splitMaster[pointi] = -2;
5110 label pointi = cPoints[i];
5112 splitMaster[pointi] = -2;
5120 label facei = mesh_.nInternalFaces();
5121 facei < mesh_.nFaces();
5125 const face&
f = mesh_.faces()[facei];
5129 splitMaster[
f[fp]] = -2;
5136 label nSplitPoints = 0;
5138 forAll(splitMaster, pointi)
5140 if (splitMaster[pointi] >= 0)
5149 forAll(splitMaster, pointi)
5151 if (splitMaster[pointi] >= 0)
5153 splitPoints[nSplitPoints++] = pointi;
5231 Pout<<
"hexRef8::consistentUnrefinement :" 5232 <<
" Determining 2:1 consistent unrefinement" <<
endl;
5238 <<
"maxSet not implemented yet." 5248 bitSet unrefinePoint(mesh_.nPoints(), pointsToUnrefine);
5255 bitSet unrefineCell(mesh_.nCells());
5257 forAll(unrefinePoint, pointi)
5259 if (unrefinePoint.test(pointi))
5261 const labelList& pCells = mesh_.pointCells(pointi);
5263 unrefineCell.
set(pCells);
5275 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
5277 label own = mesh_.faceOwner()[facei];
5278 label nei = mesh_.faceNeighbour()[facei];
5280 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5281 label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
5283 if (ownLevel < (neiLevel-1))
5290 unrefineCell.set(nei);
5301 if (!unrefineCell.test(own))
5307 unrefineCell.unset(own);
5311 else if (neiLevel < (ownLevel-1))
5315 unrefineCell.set(own);
5319 if (!unrefineCell.test(nei))
5325 unrefineCell.unset(nei);
5333 labelList neiLevel(mesh_.nBoundaryFaces());
5337 label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
5339 neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
5347 label facei = i+mesh_.nInternalFaces();
5348 label own = mesh_.faceOwner()[facei];
5349 label ownLevel = cellLevel_[own] - unrefineCell.get(own);
5351 if (ownLevel < (neiLevel[i]-1))
5355 if (!unrefineCell.test(own))
5361 unrefineCell.unset(own);
5365 else if (neiLevel[i] < (ownLevel-1))
5369 if (unrefineCell.test(own))
5375 unrefineCell.set(own);
5385 Pout<<
"hexRef8::consistentUnrefinement :" 5386 <<
" Changed " << nChanged
5387 <<
" refinement levels due to 2:1 conflicts." 5401 forAll(unrefinePoint, pointi)
5403 if (unrefinePoint.test(pointi))
5405 const labelList& pCells = mesh_.pointCells(pointi);
5409 if (!unrefineCell.test(pCells[j]))
5411 unrefinePoint.unset(pointi);
5423 forAll(unrefinePoint, pointi)
5425 if (unrefinePoint.test(pointi))
5434 forAll(unrefinePoint, pointi)
5436 if (unrefinePoint.test(pointi))
5438 newPointsToUnrefine[nSet++] = pointi;
5442 return newPointsToUnrefine;
5452 if (!history_.active())
5455 <<
"Only call if constructed with history capability" 5461 Pout<<
"hexRef8::setUnrefinement :" 5462 <<
" Checking initial mesh just to make sure" <<
endl;
5466 forAll(cellLevel_, celli)
5468 if (cellLevel_[celli] < 0)
5471 <<
"Illegal cell level " << cellLevel_[celli]
5472 <<
" for cell " << celli
5479 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5482 cellSet cSet(mesh_,
"splitPointCells", splitPointLabels.
size());
5484 forAll(splitPointLabels, i)
5486 const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
5488 cSet.insert(pCells);
5492 Pout<<
"hexRef8::setRefinement : Dumping " << pSet.
size()
5494 << cSet.size() <<
" cells for unrefinement to" <<
nl 5496 <<
" cellSet " << cSet.objectPath()
5508 forAll(splitPointLabels, i)
5512 splitFaces.insert(
pFaces);
5517 faceRemover_.compatibleRemoves
5525 if (facesToRemove.
size() != splitFaces.size())
5529 const_cast<polyMesh&
>(mesh_).setInstance
5531 mesh_.time().timeName()
5534 pointSet pSet(mesh_,
"splitPoints", splitPointLabels);
5536 faceSet fSet(mesh_,
"splitFaces", splitFaces);
5538 faceSet removeSet(mesh_,
"facesToRemove", facesToRemove);
5543 <<
"Ininitial set of split points to unrefine does not" 5544 <<
" seem to be consistent or not mid points of refined cells" 5545 <<
" splitPoints:" << splitPointLabels.
size()
5546 <<
" splitFaces:" << splitFaces.size()
5547 <<
" facesToRemove:" << facesToRemove.
size()
5556 forAll(splitPointLabels, i)
5558 label pointi = splitPointLabels[i];
5562 const labelList& pCells = mesh_.pointCells(pointi);
5565 if (pCells.
size() != 8)
5568 <<
"splitPoint " << pointi
5569 <<
" should have 8 cells using it. It has " << pCells
5578 label masterCelli =
min(pCells);
5582 label celli = pCells[j];
5584 label region = cellRegion[celli];
5589 <<
"Ininitial set of split points to unrefine does not" 5590 <<
" seem to be consistent or not mid points" 5591 <<
" of refined cells" <<
nl 5592 <<
"cell:" << celli <<
" on splitPoint " << pointi
5593 <<
" has no region to be merged into" 5597 if (masterCelli != cellRegionMaster[region])
5600 <<
"cell:" << celli <<
" on splitPoint:" << pointi
5601 <<
" in region " << region
5602 <<
" has master:" << cellRegionMaster[region]
5603 <<
" which is not the lowest numbered cell" 5604 <<
" among the pointCells:" << pCells
5613 faceRemover_.setRefinement
5624 forAll(splitPointLabels, i)
5626 label pointi = splitPointLabels[i];
5628 const labelList& pCells = mesh_.pointCells(pointi);
5630 label masterCelli =
min(pCells);
5634 cellLevel_[pCells[j]]--;
5637 history_.combineCells(masterCelli, pCells);
5641 setInstance(mesh_.facesInstance());
5651 cellLevel_.write(writeOnProc)
5652 && pointLevel_.write(writeOnProc)
5653 && level0Edge_.write(writeOnProc);
5655 if (history_.active())
5657 writeOk = writeOk && history_.write(writeOnProc);
5681 if (
exists(setsDir/
"cellLevel"))
5683 rm(setsDir/
"cellLevel");
5685 if (
exists(setsDir/
"pointLevel"))
5687 rm(setsDir/
"pointLevel");
5689 if (
exists(setsDir/
"level0Edge"))
5691 rm(setsDir/
"level0Edge");
void distribute(const mapDistributePolyMesh &)
Update local numbering for mesh redistribution.
prefixOSstream Perr
OSstream wrapped stderr (std::cerr) with parallel prefix.
void size(const label n)
Older name for setAddressableSize.
const labelList & reversePointMap() const
Reverse point map.
A class for handling file names.
readOption readOpt() const noexcept
Get the read option.
virtual const fileName & name() const
The name of the stream.
label nOldCells() const
Number of old cells.
const labelList & cellMap() const
Old cell map.
virtual label iterate(const label maxIter)
Iterate until no changes or maxIter reached.
errorManipArg< error, int > exit(error &err, const int errNo=1)
const fileName & facesInstance() const
Return the current instance directory for faces.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
void setUnrefinement(const labelList &splitPointLabels, polyTopoChange &)
Remove some refinement. Needs to be supplied output of.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void append(const T &val)
Append an element at the end of the list.
bool active() const
Is there unrefinement history?
static void removeFiles(const polyMesh &)
Helper: remove all relevant files from mesh instance.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
bool updateFace(const polyMesh &, const label thisFacei, const label neighbourCelli, const refinementData &neighbourInfo, const scalar tol, TrackingData &td)
Influence of neighbouring cell.
constexpr char nl
The newline '\n' character (0x0a)
UIndirectList< label > labelUIndList
UIndirectList of labels.
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Wave propagation of information through grid. Every iteration information goes through one layer of c...
List< cellShape > cellShapeList
List of cellShape.
void checkMesh() const
Debug: Check coupled mesh for correctness.
void storeData(const labelList &pointsToStore, const labelList &facesToStore, const labelList &cellsToStore)
Signal points/face/cells for which to store data.
Class containing mesh-to-mesh mapping information after a mesh distribution where we send parts of me...
labelList consistentSlowRefinement2(const label maxFaceDiff, const labelList &cellsToRefine, const labelList &facesToCheck) const
Like consistentSlowRefinement but uses different meshWave.
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
labelListList setRefinement(const labelList &cells, polyTopoChange &)
Insert refinement. All selected cells will be split into 8.
Holds information regarding type of cell. Used in inside/outside determination in cellClassification...
virtual const fileName & name() const override
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
List< labelList > labelListList
List of labelList.
labelList consistentUnrefinement(const labelList &pointsToUnrefine, const bool maxSet) const
Given proposed.
fileName objectPath() const
The complete path + object name.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
static scalar propagationTol() noexcept
Access to propagation tolerance.
UList< label > labelUList
A UList of labels.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
void subset(const labelList &pointMap, const labelList &faceMap, const labelList &cellMap)
Update local numbering for subsetted mesh.
fileName path() const
The complete path for the object (with instance, local,...).
#define forAll(list, i)
Loop across all elements in list.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
void setInstance(const fileName &inst)
List< face > faceList
List of faces.
label size() const noexcept
The number of elements in table.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
void clear()
Clear the list, i.e. set size to zero.
Refinement of (split) hexes using polyTopoChange.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
labelList consistentSlowRefinement(const label maxFaceDiff, const labelList &cellsToRefine, const labelList &facesToCheck, const label maxPointDiff, const labelList &pointsToCheck) const
Like consistentRefinement but slower:
Reading is optional [identical to LAZY_READ].
virtual bool write(const bool writeOnProc=true) const
Write using setting from DB.
bool exists(const fileName &name, const bool checkGzip=true, const bool followLink=true)
Does the name exist (as DIRECTORY or FILE) in the file system?
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
void distributeCellData(List< T > &values) const
Distribute list of cell data.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
Container with cells to refine. Refinement given as single direction.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
A HashTable to objects of type <T> with a labelPair key. The hashing is based on labelPair (FixedList...
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (respects is_globalIOobject trait) and check its info.
void updateMesh(const mapPolyMesh &)
Update local numbering for changed mesh.
virtual bool read()
Read object. If global number of visible cells > 0 becomes active.
const labelList & reverseCellMap() const
Reverse cell map.
int debug
Static debugging option.
Pair< label > labelPair
A pair of labels.
Type gMax(const FieldField< Field, Type > &f)
defineTypeNameAndDebug(combustionModel, 0)
const labelList & pointMap() const
Old point map.
Reduction class. If x and y are not equal assign value.
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field...
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
const labelList & visibleCells() const
Per cell in the current mesh (i.e. visible) either -1 (unrefined) or an index into splitCells...
Transfers refinement levels such that slow transition between levels is maintained. Used in FaceCellWave.
labelList consistentRefinement(const labelUList &cellLevel, const labelList &cellsToRefine, const bool maxSet) const
Given valid mesh and current cell level and proposed.
void distributePointData(List< T > &values) const
Distribute list of point data.
static void removeFiles(const polyMesh &)
Helper: remove all sets files from mesh instance.
vector point
Point is a vector.
label faceLevel(const label facei) const
Gets level such that the face has four points <= level.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
void operator()(label &x, const label y) const
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A cell is defined as a list of faces with extra functionality.
A collection of cell labels.
Direct mesh changes based on v1.3 polyTopoChange syntax.
const polyBoundaryMesh & patches
HashTable< label, labelPair, Foam::Hash< labelPair > > labelPairLookup
This is a Map of a labelPair to a label. Used for e.g. for face1, face2 to shared edge...
labelList getSplitPoints() const
Return the points at the centre of top-level split cells.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
#define DebugVar(var)
Report a variable name and value.
Mesh consisting of general polyhedral cells.
void checkRefinementLevels(const label maxPointDiff, const labelList &pointsToCheck) const
Debug: Check 2:1 consistency across faces.
List< label > labelList
A List of labels.
All refinement history. Used in unrefinement.
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, IOobject::NO_REGISTER)
bool write(const bool writeOnProc=true) const
Force writing refinement+history to polyMesh directory.
label rcIndex(const label i) const noexcept
The reverse circular index. The previous index in the list which returns to the last at the beginning...
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
ListType reorder(const labelUList &oldToNew, const ListType &input, const bool prune=false)
Reorder the elements of a list.
Defines the attributes of an object for which implicit objectRegistry management is supported...
const cellShapeList & cellShapes() const
Utility: get hexes as cell shapes.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
void setFaceInfo(const label facei, const Type &faceInfo)
Set single initial changed face.
const volScalarField & p0
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
label nOldPoints() const
Number of old points.
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
forAllConstIters(mixture.phases(), phase)
A HashTable to objects of type <T> with a label key.
static void listCombineReduce(List< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
After completion all processors have the same data.
bool rm(const fileName &file)
Remove a file (or its gz equivalent), returning true if successful.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
static constexpr const zero Zero
Global zero (0)